Information

Please refer to ‘Data Cleaning’ script prior to accessing this script.

Setup

knitr::opts_chunk$set(echo = TRUE)
require("knitr")
## Loading required package: knitr
opts_knit$set(root.dir = "~/Library/Mobile Documents/com~apple~CloudDocs/Documents/Uni/Masters/Empirical Project/Code/Empirical_Project")

# turn off scientific notation
options(scipen = 999)

## colour palette
# e7b553
# c45150
# a24b6f
# 824372
# 603863
# 382c46
# 403250

Load Libraries

library("ggplot2") # for figures
library("psych") # for Cronbach's alpha, for describe function
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library("ppcor") # for partial correlation p-values
## Loading required package: MASS
library("dplyr") # for mutate function
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library("ggpubr") # for qq-plots
library("GGally") # for scatterplot matrix 
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library("effsize") # for calculation of effect size
## 
## Attaching package: 'effsize'
## The following object is masked from 'package:psych':
## 
##     cohen.d
library("pwr") # for power calculation
library("performance") # for assessing robustness of model
library("effsize") # for eta squared
library("reshape2") # for transforming data from wide to long format
library("tidyverse") # for data cleaning 
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble  3.1.3     ✓ purrr   0.3.4
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   2.0.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x psych::%+%()    masks ggplot2::%+%()
## x psych::alpha()  masks ggplot2::alpha()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## x dplyr::select() masks MASS::select()
library("rstatix") # for ANOVA and ANCOVA
## 
## Attaching package: 'rstatix'
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
library("gridExtra") # for grid.arrange function
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
library("car") # for levene's test
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:purrr':
## 
##     some
## The following object is masked from 'package:dplyr':
## 
##     recode
## The following object is masked from 'package:psych':
## 
##     logit
library("emmeans") # to obtain estimated marginal means
## 
## Attaching package: 'emmeans'
## The following object is masked from 'package:GGally':
## 
##     pigs

Set Working Directory

# please change this to your own working directory path
setwd("~/Library/Mobile Documents/com~apple~CloudDocs/Documents/Uni/Masters/Empirical Project/Code/Empirical_Project")

Read in Data and Save Data to an Object

# please change this to however you have stored the data file 
# reading in dataframe 2, as this is the one with exclusion of n = 5
df <- read.csv(file = "data/cleaned/dataframe_2.csv", header = TRUE, na.strings = "NA")

Change Variable Classifications

# change variable classifications to meet requirements for later analyses
# ensure IVs and categorical variables are factor variables
# and DVs or continuous variables are numeric variables

# participant id and demographics
df$id <- factor(df$id)
df$age <- as.numeric(df$age)
df$sex <- factor(df$sex)
df$ethnicity <- factor(df$ethnicity)
df$sexual_orientation <- factor(df$sexual_orientation)

# fixation count DVs
df$acq_csp_fix_count <- as.numeric(df$acq_csp_fix_count)
df$acq_csm_fix_count <- as.numeric(df$acq_csm_fix_count)
df$ext_csp_fix_count <- as.numeric(df$ext_csp_fix_count)
df$ext_csm_fix_count <- as.numeric(df$ext_csm_fix_count)
df$e_ext_csp_fix_count <- as.numeric(df$e_ext_csp_fix_count)
df$l_ext_csp_fix_count <- as.numeric(df$l_ext_csp_fix_count)
df$e_ext_csm_fix_count <- as.numeric(df$e_ext_csm_fix_count)
df$l_ext_csm_fix_count <- as.numeric(df$l_ext_csm_fix_count)

# fixation duration DVs
df$acq_csp_fix_duration <- as.numeric(df$acq_csp_fix_duration)
df$acq_csm_fix_duration <- as.numeric(df$acq_csm_fix_duration)
df$ext_csp_fix_duration <- as.numeric(df$ext_csp_fix_duration)
df$ext_csm_fix_duration <- as.numeric(df$ext_csm_fix_duration)
df$e_ext_csp_fix_duration <- as.numeric(df$e_ext_csp_fix_duration)
df$l_ext_csp_fix_duration <- as.numeric(df$l_ext_csp_fix_duration)
df$e_ext_csm_fix_duration <- as.numeric(df$e_ext_csm_fix_duration)
df$l_ext_csm_fix_duration <- as.numeric(df$l_ext_csm_fix_duration)

# saccade amplitude DVs
df$acq_csp_sacc_amplitude <- as.numeric(df$acq_csp_sacc_amplitude)
df$acq_csm_sacc_amplitude <- as.numeric(df$acq_csm_sacc_amplitude)
df$ext_csp_sacc_amplitude <- as.numeric(df$ext_csp_sacc_amplitude)
df$ext_csm_sacc_amplitude <- as.numeric(df$ext_csm_sacc_amplitude)
df$e_ext_csp_sacc_amplitude <- as.numeric(df$e_ext_csp_sacc_amplitude)
df$l_ext_csp_sacc_amplitude <- as.numeric(df$l_ext_csp_sacc_amplitude)
df$e_ext_csm_sacc_amplitude <- as.numeric(df$e_ext_csm_sacc_amplitude)
df$l_ext_csm_sacc_amplitude <- as.numeric(df$l_ext_csm_sacc_amplitude)

Internal Consistency of IUS and STICSA

## IUS total
# compute & extract alpha value and save as an object
alpha_ius <- psych::alpha(df[, c("ius_1", "ius_2", "ius_3", "ius_4",
                                       "ius_5", "ius_6", "ius_7", "ius_8",
                                       "ius_9", "ius_10", "ius_11", "ius_12",
                                       "ius_13", "ius_14", "ius_15", "ius_16",
                                       "ius_17", "ius_18", "ius_19", "ius_20", 
                                       "ius_21", "ius_22", "ius_23", "ius_24",
                                       "ius_25", "ius_26", "ius_27")])$total[1]

## STICSA total
# compute & extract alpha value and save as an object
alpha_sticsa <- psych::alpha(df[, c("sticsa_1", "sticsa_2", "sticsa_3", "sticsa_4",
                                       "sticsa_5", "sticsa_6", "sticsa_7", "sticsa_8",
                                       "sticsa_9", "sticsa_10", "sticsa_11", "sticsa_12",
                                       "sticsa_13", "sticsa_14", "sticsa_15", "sticsa_16",
                                       "sticsa_17", "sticsa_18", "sticsa_19", "sticsa_20", 
                                       "sticsa_21")])$total[1]

# create table of both Crobach's alpha values
cronbachs_alpha_questionnaires <- rbind(alpha_ius, alpha_sticsa)

# clean up row and column names for easier interpretation
rownames(cronbachs_alpha_questionnaires) <- c("IUS-27", "STICSA")
colnames(cronbachs_alpha_questionnaires) <- "Cronbach's Alpha"

# obtain Cronbach's alpha table
cronbachs_alpha_questionnaires
##        Cronbach's Alpha
## IUS-27        0.9496736
## STICSA        0.8766597

Compute Questionnaire Totals

#### IUS total
# all items, no reverse scoring
df$ius_total <- as.numeric(df$ius_1 + df$ius_2 + df$ius_3 + df$ius_4 + df$ius_5 + df$ius_6 + 
  df$ius_7 + df$ius_8 + df$ius_9 + df$ius_10 + df$ius_11 + df$ius_12 + df$ius_13 + 
  df$ius_14 + df$ius_15 + df$ius_16 + df$ius_17 + df$ius_18 + df$ius_19 + df$ius_20 + 
  df$ius_21 + df$ius_22 + df$ius_23 + df$ius_24 + df$ius_25 + df$ius_26 + df$ius_27)

#### STICSA total
# all items, no reverse scoring
df$sticsa_total <- as.numeric(df$sticsa_1 + df$sticsa_2 + df$sticsa_3 +
                                df$sticsa_4 + df$sticsa_5 + df$sticsa_6 +
                                df$sticsa_7 + df$sticsa_8 + df$sticsa_9 +
                                df$sticsa_10 + df$sticsa_11 + df$sticsa_12 +
                                df$sticsa_13 + df$sticsa_14 + df$sticsa_15 +
                                df$sticsa_16 + df$sticsa_17 + df$sticsa_18 +
                                df$sticsa_19 + df$sticsa_20 + df$sticsa_21)

Create High / Low IU Classifications

# compute variable classifying participants as high/ low IU on basis of median split, 
# and store as factor
df$iu_group <- factor(ifelse(df$ius_total >= 65, 1, -1))
# high IU = 1
# low IU = -1

Check Distribution and Range to Identify Extreme Scores and Potential Data Errors in Questionnaires

For IUS 27 Total in Both Groups

# possible total scores for the IUS range from 27-135 

########################## check distributions 
hist_ius_total <- df %>%
  ggplot(aes(ius_total, fill = iu_group)) +
  geom_histogram(binwidth = 5, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(20, 140, 10)) +
  labs(x = "IUS-27 Total", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram of IUS-27 Scores") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_ius_total

# save plot to file
ggsave(filename = "graphs/histograms/hist_ius_total.png", 
       plot = hist_ius_total, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

########################## check ranges
range_ius_total <- by(df$ius_total, df$iu_group, range)
range_ius_total 
## df$iu_group: -1
## [1] 32 64
## ------------------------------------------------------------ 
## df$iu_group: 1
## [1]  65 125
# for high IU: 65-125
# for low IU: 32-64
##### overall: all scores are in range of possible scores, no errors apparent

For STICSA Total in Both Groups

# possible total scores for the STICSA range from 21-84

########################## check distributions 
hist_sticsa_total <- df %>%
  ggplot(aes(sticsa_total, fill = iu_group)) +
  geom_histogram(binwidth = 5, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(20, 90, 10)) + 
  labs(x = "STICSA Total", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram of STICSA Scores") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_sticsa_total

# save plot to file
ggsave(filename = "graphs/histograms/hist_sticsa_total.png", 
       plot = hist_sticsa_total, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

########################## check ranges
range_sticsa_total <- by(df$sticsa_total, df$iu_group, range)
range_sticsa_total 
## df$iu_group: -1
## [1] 22 57
## ------------------------------------------------------------ 
## df$iu_group: 1
## [1] 30 69
# for high IU: 30-69
# for low IU: 22-57
##### overall: all scores are in range of possible scores, no errors apparent

Compute Demographics

#### for age
# for all participants
all_age_table <- 
  describe(df[, "age"])

# for high IU
high_iu_age_table <- 
  describe(df[df$iu_group =="1", "age"])

# for low IU
low_iu_age_table <- 
  describe(df[df$iu_group =="-1", "age"])

# combine in a table
age_table <- rbind(all_age_table, high_iu_age_table, low_iu_age_table)

# re-name rows for easier interpretation
rownames(age_table) <- c("Age (All Participants","Age (High IU Group)", 
                         "Age (Low IU Group)")

### for sex 
sex_table <- xtabs(~ iu_group + sex, data = df)
sex_table <- prop.table(sex_table) %>%
  round(digits = 4) * 100
rownames(sex_table) <- c("Low IU", "High IU")
sex_table
##          sex
## iu_group  Female  Male
##   Low IU   26.28 24.82
##   High IU  34.31 14.60
### for sexual orientation
sexual_orientation_table <- xtabs(~ iu_group + sexual_orientation, data = df)
sexual_orientation_table <- prop.table(sexual_orientation_table) %>%
  round(digits = 4) * 100
rownames(sexual_orientation_table) <- c("Low IU", "High IU")
sexual_orientation_table
##          sexual_orientation
## iu_group  Heterosexual Sexual Minority
##   Low IU         42.15            7.44
##   High IU        42.98            7.44
### for ethnicity
ethnicity_table <- xtabs(~ iu_group + ethnicity, data = df)
ethnicity_table <- prop.table(ethnicity_table) %>%
  round(digits = 4) * 100
rownames(ethnicity_table) <- c("Low IU", "High IU")
ethnicity_table
##          ethnicity
## iu_group  Asian Black Middle Eastern/ Arab Mixed White
##   Low IU   7.26  1.61                 2.42  0.81 37.90
##   High IU 16.13  0.00                 0.81  0.81 32.26
#### write each to csv 
# age
write.csv(age_table, file = "tables/demographics/age_table.csv",
          row.names = TRUE)

# ethnicity
write.csv(ethnicity_table, file = "tables/demographics/ethnicity_table.csv",
          row.names = TRUE)

# sex
write.csv(sex_table, file = "tables/demographics/sex_table.csv", 
          row.names = TRUE)

# sexual orientation
write.csv(sexual_orientation_table, file = "tables/demographics/sexual_orientation_table.csv",
          row.names = TRUE)

Check for Difference in Questionnaire Totals Between Groups

Check for Difference in IUS-27 Totals Between Groups

# t-test to check for intergroup differences in IU total score - there should 
# be a significant difference. 

# plot data for both groups using QQ plot
qqplot_ius_total <- ggqqplot(df, x = "ius_total",
         color = "iu_group",
         palette = c("#c45150", "#824372"),
         title = "Q-Q Plot IUS Total") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) 
# inspect the QQ plots
qqplot_ius_total 

# save plot to file
ggsave(filename = "graphs/qqplots/qqplot_ius_total.png", 
       plot = qqplot_ius_total, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

# check significance of data for both groups using Shapiro-Wilk Test
shapiro_ius_total <- by(df$ius_total, df$iu_group, shapiro.test)
shapiro_ius_total
## df$iu_group: -1
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.9603, p-value = 0.02444
## 
## ------------------------------------------------------------ 
## df$iu_group: 1
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.88758, p-value = 0.00001673
# high IU: p-value < .05, data violate assumption of normality
# low IU: p-value < .05, data violate assumption of normality

## check assumption of homogeneity of variances using Bartlett Test ##
bartlett_ius_total <- bartlett.test(ius_total ~ iu_group, data = df) 
bartlett_ius_total
## 
##  Bartlett test of homogeneity of variances
## 
## data:  ius_total by iu_group
## Bartlett's K-squared = 19.86, df = 1, p-value = 0.000008334
# p-value < .05, data violate assumption of equal variances

## compute independent samples t.test ##
# as data violate assumption of normality and assumption of equal variances,
# use non-parametric Mann Whitney U

# compute t.test and assign values to an object
ius_total_groupdiff <- wilcox.test(ius_total ~ iu_group, data = df, paired = FALSE)
# obtain t.test values 
ius_total_groupdiff
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  ius_total by iu_group
## W = 0, p-value < 0.00000000000000022
## alternative hypothesis: true location shift is not equal to 0
# p-value < .05, there is a statistical difference in IUS 27 total between groups

Check for Difference in STICSA Totals Between Groups

# t-test to check for intergroup differences in STICSA total score - there should 
# be a significant difference. 

# plot data for both groups using QQ plot
qqplot_sticsa_total <- ggqqplot(df, x = "sticsa_total",
         color = "iu_group",
         palette = c("#c45150", "#824372"),
         title = "Q-Q Plot STICSA Total") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) 
# inspect the QQ plots
qqplot_sticsa_total 

# save plot to file
ggsave(filename = "graphs/qqplots/qqplot_sticsa_total.png", 
       plot = qqplot_sticsa_total, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

# check significance of data for both groups using Shapiro-Wilk Test
shapiro_sticsa_total <- by(df$sticsa_total, df$iu_group, shapiro.test)
shapiro_sticsa_total
## df$iu_group: -1
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.95514, p-value = 0.01269
## 
## ------------------------------------------------------------ 
## df$iu_group: 1
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.96452, p-value = 0.05002
# high iu: p-value > .05, data meet assumption of normality
# low iu: p-value < .05, data violate assumption of normality

## check assumption of homogeneity of variances using Bartlett Test ##
bartlett_sticsa_total <- bartlett.test(sticsa_total ~ iu_group, data = df) 
bartlett_sticsa_total
## 
##  Bartlett test of homogeneity of variances
## 
## data:  sticsa_total by iu_group
## Bartlett's K-squared = 3.8861, df = 1, p-value = 0.04869
# p-value < .05, data violate assumption of equal variances

## compute independent samples t.test ##
# as data violate assumption of normality and assumption of equal variances,
# use non-parametric Mann Whitney U

# compute t.test and assign values to an object
sticsa_total_groupdiff <- wilcox.test(sticsa_total ~ iu_group, data = df, paired = FALSE)
# obtain t.test values 
sticsa_total_groupdiff
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  sticsa_total by iu_group
## W = 1032, p-value = 0.000000005691
## alternative hypothesis: true location shift is not equal to 0
# p-value < .05, there is a statistical difference in STICSA total between groups

Check for Difference in Demographics Between Groups

Check for Difference in Age Between Groups

# t-test to check for intergroup differences in age

# plot data for both groups using QQ plot
qqplot_age <- ggqqplot(df, x = "age",
         color = "iu_group",
         palette = c("#c45150", "#824372"),
         title = "Q-Q Plot Age") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) 
# inspect the QQ plots
qqplot_age 
## Warning: Removed 1 rows containing non-finite values (stat_qq).
## Warning: Removed 1 rows containing non-finite values (stat_qq_line).

## Warning: Removed 1 rows containing non-finite values (stat_qq_line).

# save plot to file
ggsave(filename = "graphs/qqplots/qqplot_age.png", 
       plot = qqplot_age, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")
## Warning: Removed 1 rows containing non-finite values (stat_qq).

## Warning: Removed 1 rows containing non-finite values (stat_qq_line).

## Warning: Removed 1 rows containing non-finite values (stat_qq_line).
# check significance of data for both groups using Shapiro-Wilk Test
shapiro_age <- by(df$age, df$iu_group, shapiro.test)
shapiro_age
## df$iu_group: -1
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.95698, p-value = 0.016
## 
## ------------------------------------------------------------ 
## df$iu_group: 1
## 
##  Shapiro-Wilk normality test
## 
## data:  dd[x, ]
## W = 0.88408, p-value = 0.00001422
# high IU: p-value < .05, data violate assumption of normality
# low IU: p-value < .05, data violate assumption of normality

## check assumption of homogeneity of variances using Bartlett Test ##
bartlett_age <- bartlett.test(age ~ iu_group, data = df) 
bartlett_age
## 
##  Bartlett test of homogeneity of variances
## 
## data:  age by iu_group
## Bartlett's K-squared = 0.27665, df = 1, p-value = 0.5989
# p-value > .05, data meet assumption of equal variances

## compute independent samples t.test ##
# as data violate assumption of normality,
# use non-parametric Mann Whitney U

# compute t.test and assign values to an object
age_groupdiff <- wilcox.test(age ~ iu_group, data = df, paired = FALSE)

# obtain t.test values 
age_groupdiff
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  age by iu_group
## W = 2585.5, p-value = 0.3773
## alternative hypothesis: true location shift is not equal to 0
# p-value > .05, there is no statistical difference in age between groups

Check for Difference in Sex Between Groups

# compute chi-square of cross-tabulation and save as object 
chi_sex <- chisq.test(table(df$iu_group, df$sex))

# check assumption of chi-square
chi_sex_expected <- chi_sex$expected 
chi_sex_expected
##     
##        Female     Male
##   -1 42.40876 27.59124
##   1  40.59124 26.40876
# no cells less than 5, meets assumptions 

# obtain statistic, df and p-value
chi_sex
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(df$iu_group, df$sex)
## X-squared = 4.2708, df = 1, p-value = 0.03877
# p-value < .05, there appears to be a statistical difference in sex between groups

# therefore, obtain observed values
chi_sex_observed <- chi_sex$observed
chi_sex_observed
##     
##      Female Male
##   -1     36   34
##   1      47   20

Check for Difference in Ethnicity Between Groups

# compute chi-square of cross-tabulation and save as object 
chi_ethnicity <- chisq.test(table(df$iu_group, df$ethnicity)) 
## Warning in chisq.test(table(df$iu_group, df$ethnicity)): Chi-squared
## approximation may be incorrect
# check assumption of chi-square
chi_ethnicity$expected 
##     
##      Asian Black Middle Eastern/ Arab Mixed White
##   -1  14.5     1                    2     1  43.5
##   1   14.5     1                    2     1  43.5
# multiple cells with values less than 5, does not meet assumptions
# and therefore requires Fisher's Exact Test

# obtain statistic and df
chi_ethnicity 
## 
##  Pearson's Chi-squared test
## 
## data:  table(df$iu_group, df$ethnicity)
## X-squared = 7.7356, df = 4, p-value = 0.1018
# obtain corrected p-value 
chi_ethnicity_pval <- fisher.test(df$iu_group, df$ethnicity) 
chi_ethnicity_pval
## 
##  Fisher's Exact Test for Count Data
## 
## data:  df$iu_group and df$ethnicity
## p-value = 0.05899
## alternative hypothesis: two.sided
# p-value > .05, no evidence of statistical difference in ethnicity between groups

Check for Difference in Sexual Orientation Between Groups

# compute chi-square of cross-tabulation and save as object 
chi_sexual_orientation <- chisq.test(table(df$iu_group, df$sexual_orientation)) 

# check assumption of chi-square
chi_sexual_orientation$expected 
##     
##      Heterosexual Sexual Minority
##   -1     51.07438         8.92562
##   1      51.92562         9.07438
# no cells with values less than 5, meets assumptions

# obtain statistic and df
chi_sexual_orientation 
## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(df$iu_group, df$sexual_orientation)
## X-squared = 0, df = 1, p-value = 1
# p-value > .05, no evidence of statistical difference in sexual orientation between groups

Distribution Checks of Eye-Movement Variables

Fixation Count

Acquisition CS+

hist_acq_csp_fix_count <- df %>%
  ggplot(aes(acq_csp_fix_count, fill = iu_group)) +
  geom_histogram(binwidth = 1, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 30, 5)) +
  labs(x = "Fixation Count", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Acquisition CS+ Fixation Count") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_acq_csp_fix_count

# save plot to file
ggsave(filename = "graphs/histograms/hist_acq_csp_fix_count.png", 
       plot = hist_acq_csp_fix_count, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Acquisition CS-

hist_acq_csm_fix_count <- df %>%
  ggplot(aes(acq_csm_fix_count, fill = iu_group)) +
  geom_histogram(binwidth = 1, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 30, 5)) +
  labs(x = "Fixation Count", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Acquisition CS- Fixation Count") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_acq_csm_fix_count

# save plot to file
ggsave(filename = "graphs/histograms/hist_acq_csm_fix_count.png", 
       plot = hist_acq_csm_fix_count, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

# combine acquisition fixation count graphs
hists_acq_fix_count <- 
  grid.arrange(hist_acq_csp_fix_count, hist_acq_csm_fix_count,
             ncol =2)

# save plot to file
ggsave(filename = "graphs/histograms/hists_acq_fix_count.png", 
       plot = hists_acq_fix_count, 
       width = 30, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Early Extinction CS+

hist_e_ext_csp_fix_count <- df %>%
  ggplot(aes(e_ext_csp_fix_count, fill = iu_group)) +
  geom_histogram(binwidth = 1, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 30, 5)) +
  labs(x = "Fixation Count", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Early Extinction CS+ Fixation Count") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_e_ext_csp_fix_count

# save plot to file
ggsave(filename = "graphs/histograms/hist_e_ext_csp_fix_count.png", 
       plot = hist_e_ext_csp_fix_count, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Early Extinction CS-

hist_e_ext_csm_fix_count <- df %>%
  ggplot(aes(e_ext_csm_fix_count, fill = iu_group)) +
  geom_histogram(binwidth = 1, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 30, 5)) +
  labs(x = "Fixation Count", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Early Extinction CS- Fixation Count") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_e_ext_csm_fix_count

# save plot to file
ggsave(filename = "graphs/histograms/hist_e_ext_csm_fix_count.png", 
       plot = hist_e_ext_csm_fix_count, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Late Extinction CS+

hist_l_ext_csp_fix_count <- df %>%
  ggplot(aes(l_ext_csp_fix_count, fill = iu_group)) +
  geom_histogram(binwidth = 1, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 30, 5)) +
  labs(x = "Fixation Count", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Late Extinction CS+ Fixation Count") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_l_ext_csp_fix_count

# save plot to file
ggsave(filename = "graphs/histograms/hist_l_ext_csp_fix_count.png", 
       plot = hist_l_ext_csp_fix_count, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Late Extinction CS-

hist_l_ext_csm_fix_count <- df %>%
  ggplot(aes(l_ext_csm_fix_count, fill = iu_group)) +
  geom_histogram(binwidth = 1, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 30, 5)) +
  labs(x = "Fixation Count", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Late Extinction CS- Fixation Count") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_l_ext_csm_fix_count

# save plot to file
ggsave(filename = "graphs/histograms/hist_l_ext_csm_fix_count.png", 
       plot = hist_l_ext_csm_fix_count, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

# combine extinction fixation count graphs
hists_ext_fix_count <- 
  grid.arrange(hist_e_ext_csp_fix_count, hist_e_ext_csm_fix_count,
               hist_l_ext_csp_fix_count, hist_l_ext_csm_fix_count,
             ncol =2)

# save plot to file
ggsave(filename = "graphs/histograms/hists_ext_fix_count.png", 
       plot = hists_ext_fix_count, 
       width = 30, 
       height = 20, 
       dpi = 300, 
       units = "cm")

Fixation Duration

Acquisition CS+

hist_acq_csp_fix_duration <- df %>%
  ggplot(aes(acq_csp_fix_duration, fill = iu_group)) +
  geom_histogram(binwidth = 220, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 6000, 1000)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Acquisition CS+ Fixation Duration") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_acq_csp_fix_duration

# save plot to file
ggsave(filename = "graphs/histograms/hist_acq_csp_fix_duration.png", 
       plot = hist_acq_csp_fix_duration, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Acquisition CS-

hist_acq_csm_fix_duration <- df %>%
  ggplot(aes(acq_csm_fix_duration, fill = iu_group)) +
  geom_histogram(binwidth = 220, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 6000, 1000)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Acquisition CS- Fixation Duration") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_acq_csm_fix_duration

# save plot to file
ggsave(filename = "graphs/histograms/hist_acq_csm_fix_duration.png", 
       plot = hist_acq_csm_fix_duration, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

# combine acquisition fixation duration graphs
hists_acq_fix_duration <- 
  grid.arrange(hist_acq_csp_fix_duration, hist_acq_csm_fix_duration,
             ncol =2)

# save plot to file
ggsave(filename = "graphs/histograms/hists_acq_fix_duration.png", 
       plot = hists_acq_fix_duration, 
       width = 30, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Early Extinction CS+

hist_e_ext_csp_fix_duration <- df %>%
  ggplot(aes(e_ext_csp_fix_duration, fill = iu_group)) +
  geom_histogram(binwidth = 220, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 6000, 1000)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Early Extinction CS+ Fixation Duration") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_e_ext_csp_fix_duration

# save plot to file
ggsave(filename = "graphs/histograms/hist_e_ext_csp_fix_duration.png", 
       plot = hist_e_ext_csp_fix_duration, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Early Extinction CS-

hist_e_ext_csm_fix_duration <- df %>%
  ggplot(aes(e_ext_csm_fix_duration, fill = iu_group)) +
  geom_histogram(binwidth = 220, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 6000, 1000)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Early Extinction CS- Fixation Duration") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_e_ext_csm_fix_duration

# save plot to file
ggsave(filename = "graphs/histograms/hist_e_ext_csm_fix_duration.png", 
       plot = hist_e_ext_csm_fix_duration, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Late Extinction CS+

hist_l_ext_csp_fix_duration <- df %>%
  ggplot(aes(l_ext_csp_fix_duration, fill = iu_group)) +
  geom_histogram(binwidth = 220, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 6000, 1000)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Late Extinction CS+ Fixation Duration") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_l_ext_csp_fix_duration

# save plot to file
ggsave(filename = "graphs/histograms/hist_l_ext_csp_fix_duration.png", 
       plot = hist_l_ext_csp_fix_duration, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Late Extinction CS-

hist_l_ext_csm_fix_duration <- df %>%
  ggplot(aes(l_ext_csm_fix_duration, fill = iu_group)) +
  geom_histogram(binwidth = 220, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 6000, 1000)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Late Extinction CS- Fixation Duration") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_l_ext_csm_fix_duration

# save plot to file
ggsave(filename = "graphs/histograms/hist_l_ext_csm_fix_duration.png", 
       plot = hist_l_ext_csm_fix_duration, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

# combine extinction fixation duration graphs
hists_ext_fix_duration <- 
  grid.arrange(hist_e_ext_csp_fix_duration, hist_e_ext_csm_fix_duration,
               hist_l_ext_csp_fix_duration, hist_l_ext_csm_fix_duration,
             ncol =2)

# save plot to file
ggsave(filename = "graphs/histograms/hists_ext_fix_duration.png", 
       plot = hists_ext_fix_duration, 
       width = 30, 
       height = 20, 
       dpi = 300, 
       units = "cm")

Saccade Amplitude

Acquisition CS+

hist_acq_csp_sacc_amplitude <- df %>%
  ggplot(aes(acq_csp_sacc_amplitude, fill = iu_group)) +
  geom_histogram(binwidth = .5, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  labs(x = "Saccade Amplitude (deg/ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Acquisition CS+ Saccade Amplitude") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_acq_csp_sacc_amplitude

ggsave(filename = "graphs/histograms/hist_acq_csp_sacc_amplitude.png", 
       plot = hist_acq_csp_sacc_amplitude, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Acquisition CS-

hist_acq_csm_sacc_amplitude <- df %>%
  ggplot(aes(acq_csm_sacc_amplitude, fill = iu_group)) +
  geom_histogram(binwidth = .5, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 10, 2)) +
  labs(x = "Saccade Amplitude (deg/ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Acquisition CS- Saccade Amplitude") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_acq_csm_sacc_amplitude
## Warning: Removed 2 rows containing non-finite values (stat_bin).

ggsave(filename = "graphs/histograms/hist_acq_csm_sacc_amplitude.png", 
       plot = hist_acq_csm_sacc_amplitude, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")
## Warning: Removed 2 rows containing non-finite values (stat_bin).
# combine acquisition saccade amplitude graphs
hists_acq_sacc_amplitude <- 
  grid.arrange(hist_acq_csp_sacc_amplitude, hist_acq_csm_sacc_amplitude,
             ncol =2)
## Warning: Removed 2 rows containing non-finite values (stat_bin).

# save plot to file
ggsave(filename = "graphs/histograms/hists_acq_sacc_amplitude.png", 
       plot = hists_acq_sacc_amplitude, 
       width = 30, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Early Extinction CS+

hist_e_ext_csp_sacc_amplitude <- df %>%
  ggplot(aes(e_ext_csp_sacc_amplitude, fill = iu_group)) +
  geom_histogram(binwidth = .5, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 14, 2)) +
  labs(x = "Saccade Amplitude (deg/ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Early Extinction CS+ Saccade Amplitude") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_e_ext_csp_sacc_amplitude
## Warning: Removed 1 rows containing non-finite values (stat_bin).

# save plot to file
ggsave(filename = "graphs/histograms/hist_e_ext_csp_sacc_amplitude.png", 
       plot = hist_e_ext_csp_sacc_amplitude, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")
## Warning: Removed 1 rows containing non-finite values (stat_bin).

Early Extinction CS-

hist_e_ext_csm_sacc_amplitude <- df %>%
  ggplot(aes(e_ext_csm_sacc_amplitude, fill = iu_group)) +
  geom_histogram(binwidth = .5, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 14, 2)) +
  labs(x = "Saccade Amplitude (deg/ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Early Extinction CS- Saccade Amplitude") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_e_ext_csm_sacc_amplitude
## Warning: Removed 1 rows containing non-finite values (stat_bin).

# save plot to file
ggsave(filename = "graphs/histograms/hist_e_ext_csm_sacc_amplitude.png", 
       plot = hist_e_ext_csm_sacc_amplitude, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")
## Warning: Removed 1 rows containing non-finite values (stat_bin).

Late Extinction CS+

hist_l_ext_csp_sacc_amplitude <- df %>%
  ggplot(aes(l_ext_csp_sacc_amplitude, fill = iu_group)) +
  geom_histogram(binwidth = .5, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 14, 2)) +
  labs(x = "Saccade Amplitude (deg/ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Late Extinction CS+ Saccade Amplitude") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_l_ext_csp_sacc_amplitude
## Warning: Removed 1 rows containing non-finite values (stat_bin).

# save plot to file
ggsave(filename = "graphs/histograms/hist_l_ext_csp_sacc_amplitude.png", 
       plot = hist_l_ext_csp_sacc_amplitude, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")
## Warning: Removed 1 rows containing non-finite values (stat_bin).

Late Extinction CS-

hist_l_ext_csm_sacc_amplitude <- df %>%
  ggplot(aes(l_ext_csm_sacc_amplitude, fill = iu_group)) +
  geom_histogram(binwidth = .5, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 14, 2)) +
  labs(x = "Saccade Amplitude (deg/ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Late Extinction CS- Saccade Amplitude") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_l_ext_csm_sacc_amplitude

# save plot to file
ggsave(filename = "graphs/histograms/hist_l_ext_csm_sacc_amplitude.png", 
       plot = hist_l_ext_csm_sacc_amplitude, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

# combine extinction saccade amplitude graphs
hists_ext_sacc_amplitude <- 
  grid.arrange(hist_e_ext_csp_sacc_amplitude, hist_e_ext_csm_sacc_amplitude,
               hist_l_ext_csp_sacc_amplitude, hist_l_ext_csm_sacc_amplitude,
             ncol =2)
## Warning: Removed 1 rows containing non-finite values (stat_bin).

## Warning: Removed 1 rows containing non-finite values (stat_bin).

## Warning: Removed 1 rows containing non-finite values (stat_bin).

# save plot to file
ggsave(filename = "graphs/histograms/hists_ext_sacc_amplitude.png", 
       plot = hists_ext_sacc_amplitude, 
       width = 30, 
       height = 20, 
       dpi = 300, 
       units = "cm")

Descriptives

Questionnaire Variables

# for all participants
descriptives_all_questionnaires <- 
  describe(df[, c("ius_total", "sticsa_total")], na.rm = TRUE)

# for high IU group
descriptives_high_iu_questionnaires <- 
  describe(df[df$iu_group == "1", c("ius_total", "sticsa_total")], na.rm = TRUE)

# for low IU group
descriptives_low_iu_questionnaires <- 
  describe(df[df$iu_group == "-1", c("ius_total", "sticsa_total")], na.rm = TRUE)

# combine all into table
descriptives_questionnaires_table <- round(rbind(descriptives_all_questionnaires, 
                                           descriptives_high_iu_questionnaires,
                                           descriptives_low_iu_questionnaires), 2)

# rename rows for easier interpretation
rownames(descriptives_questionnaires_table) <- c("IUS 27 (All Participants)",
                                                 "STICSA Total (All Participants)",
                                                 "IUS 27 (High IU Group)",
                                                 "STICSA Total (High IU Group)",
                                                 "IUS 27 (Low IU Group)",
                                                 "STICSA Total (Low IU Group)")

descriptives_questionnaires_table
##                                 vars   n  mean    sd median trimmed   mad min
## IUS 27 (All Participants)          1 139 65.82 20.39   63.0   64.27 20.76  32
## STICSA Total (All Participants)    2 139 40.54  9.54   39.0   39.93 10.38  22
## IUS 27 (High IU Group)             1  68 82.65 14.77   78.0   80.79 11.86  65
## STICSA Total (High IU Group)       2  68 45.29  9.30   45.5   44.77  9.64  30
## IUS 27 (Low IU Group)              1  71 49.70  8.51   51.0   49.96 10.38  32
## STICSA Total (Low IU Group)        2  71 35.99  7.32   35.0   35.35  5.93  22
##                                 max range  skew kurtosis   se
## IUS 27 (All Participants)       125    93  0.64     0.00 1.73
## STICSA Total (All Participants)  69    47  0.65     0.06 0.81
## IUS 27 (High IU Group)          125    60  1.11     0.53 1.79
## STICSA Total (High IU Group)     69    39  0.47    -0.19 1.13
## IUS 27 (Low IU Group)            64    32 -0.23    -1.05 1.01
## STICSA Total (Low IU Group)      57    35  0.76     0.30 0.87
# write to csv
write.csv(descriptives_questionnaires_table, file = "tables/descriptives/descriptives_questionnaires_table.csv",
          row.names = TRUE)

Eye Movement Variables

Fixation Count

# for all participants
descriptives_all_fix_count <- 
  describe(df[, c("acq_csp_fix_count","acq_csm_fix_count", 
                  "e_ext_csp_fix_count", "e_ext_csm_fix_count",
                  "l_ext_csp_fix_count", "l_ext_csm_fix_count")], 
           na.rm = TRUE)
  
# for high IU group
descriptives_high_iu_fix_count <- 
  describe(df[df$iu_group == "1", c("acq_csp_fix_count","acq_csm_fix_count", 
                                    "e_ext_csp_fix_count", "e_ext_csm_fix_count",
                                    "l_ext_csp_fix_count", "l_ext_csm_fix_count")], 
           na.rm = TRUE)
  
# for low IU group
descriptives_low_iu_fix_count <- 
  describe(df[df$iu_group == "-1", c("acq_csp_fix_count","acq_csm_fix_count", 
                                     "e_ext_csp_fix_count", "e_ext_csm_fix_count",
                                     "l_ext_csp_fix_count", "l_ext_csm_fix_count")], 
           na.rm = TRUE)

# combine all into table
descriptives_fix_count_table <- round(rbind(descriptives_all_fix_count, 
                                            descriptives_high_iu_fix_count,
                                            descriptives_low_iu_fix_count), 2)

# rename rows for easier interpretation
rownames(descriptives_fix_count_table) <- c("Acquisition CS+ Fix Count (All Participants)", 
                                            "Acquisition CS- Fix Count (All Participants)", 
                                            "Early Extinction CS+ Fix Count (All Participants)", 
                                            "Early Extinction CS- Fix Count (All Participants)",
                                            "Late Extinction CS+ Fix Count (All Participants)", 
                                            "Late Extinction CS- Fix Count (All Participants)",
                                            "Acquisition CS+ Fix Count (High IU Group)", 
                                            "Acquisition CS- Fix Count (High IU Group)", 
                                            "Early Extinction CS+ Fix Count (High IU Group)", 
                                            "Early Extinction CS- Fix Count (High IU Group)",
                                            "Late Extinction CS+ Fix Count (High IU Group)", 
                                            "Late Extinction CS- Fix Count (High IU Group)",
                                            "Acquisition CS+ Fix Count (Low IU Group)", 
                                            "Acquisition CS- Fix Count (Low IU Group)", 
                                            "Early Extinction CS+ Fix Count (Low IU Group)", 
                                            "Early Extinction CS- Fix Count (Low IU Group)",
                                            "Late Extinction CS+ Fix Count (Low IU Group)", 
                                            "Late Extinction CS- Fix Count (Low IU Group)")

descriptives_fix_count_table
##                                                   vars   n mean   sd median
## Acquisition CS+ Fix Count (All Participants)         1 139 6.90 3.65   6.67
## Acquisition CS- Fix Count (All Participants)         2 139 7.31 3.25   6.75
## Early Extinction CS+ Fix Count (All Participants)    3 139 7.16 3.70   6.50
## Early Extinction CS- Fix Count (All Participants)    4 139 7.40 3.53   6.75
## Late Extinction CS+ Fix Count (All Participants)     5 139 7.55 3.49   7.25
## Late Extinction CS- Fix Count (All Participants)     6 139 7.86 3.52   7.75
## Acquisition CS+ Fix Count (High IU Group)            1  68 7.51 3.84   7.08
## Acquisition CS- Fix Count (High IU Group)            2  68 7.97 3.07   7.79
## Early Extinction CS+ Fix Count (High IU Group)       3  68 7.54 3.26   6.75
## Early Extinction CS- Fix Count (High IU Group)       4  68 8.14 3.26   7.88
## Late Extinction CS+ Fix Count (High IU Group)        5  68 8.41 3.63   7.75
## Late Extinction CS- Fix Count (High IU Group)        6  68 8.89 3.33   8.75
## Acquisition CS+ Fix Count (Low IU Group)             1  71 6.33 3.38   5.50
## Acquisition CS- Fix Count (Low IU Group)             2  71 6.67 3.31   5.92
## Early Extinction CS+ Fix Count (Low IU Group)        3  71 6.80 4.06   6.00
## Early Extinction CS- Fix Count (Low IU Group)        4  71 6.70 3.66   5.75
## Late Extinction CS+ Fix Count (Low IU Group)         5  71 6.72 3.15   6.50
## Late Extinction CS- Fix Count (Low IU Group)         6  71 6.87 3.43   6.50
##                                                   trimmed  mad  min   max range
## Acquisition CS+ Fix Count (All Participants)         6.57 3.71 1.50 23.17 21.67
## Acquisition CS- Fix Count (All Participants)         7.06 3.71 1.92 18.33 16.42
## Early Extinction CS+ Fix Count (All Participants)    6.75 3.71 1.50 20.50 19.00
## Early Extinction CS- Fix Count (All Participants)    7.14 3.71 1.50 21.50 20.00
## Late Extinction CS+ Fix Count (All Participants)     7.33 3.34 1.00 22.00 21.00
## Late Extinction CS- Fix Count (All Participants)     7.65 3.34 1.50 20.00 18.50
## Acquisition CS+ Fix Count (High IU Group)            7.14 3.46 2.00 23.17 21.17
## Acquisition CS- Fix Count (High IU Group)            7.78 3.21 2.92 18.33 15.42
## Early Extinction CS+ Fix Count (High IU Group)       7.26 2.97 2.25 19.25 17.00
## Early Extinction CS- Fix Count (High IU Group)       8.00 3.34 2.00 16.50 14.50
## Late Extinction CS+ Fix Count (High IU Group)        8.11 3.71 1.50 22.00 20.50
## Late Extinction CS- Fix Count (High IU Group)        8.61 2.41 3.25 20.00 16.75
## Acquisition CS+ Fix Count (Low IU Group)             6.02 3.71 1.50 15.67 14.17
## Acquisition CS- Fix Count (Low IU Group)             6.35 3.21 1.92 15.50 13.58
## Early Extinction CS+ Fix Count (Low IU Group)        6.22 3.71 1.50 20.50 19.00
## Early Extinction CS- Fix Count (Low IU Group)        6.31 3.34 1.50 21.50 20.00
## Late Extinction CS+ Fix Count (Low IU Group)         6.57 3.34 1.00 16.00 15.00
## Late Extinction CS- Fix Count (Low IU Group)         6.68 3.71 1.50 17.75 16.25
##                                                   skew kurtosis   se
## Acquisition CS+ Fix Count (All Participants)      1.18     2.52 0.31
## Acquisition CS- Fix Count (All Participants)      0.69     0.19 0.28
## Early Extinction CS+ Fix Count (All Participants) 1.18     1.74 0.31
## Early Extinction CS- Fix Count (All Participants) 0.82     0.96 0.30
## Late Extinction CS+ Fix Count (All Participants)  0.83     1.34 0.30
## Late Extinction CS- Fix Count (All Participants)  0.68     0.81 0.30
## Acquisition CS+ Fix Count (High IU Group)         1.43     3.48 0.47
## Acquisition CS- Fix Count (High IU Group)         0.74     0.76 0.37
## Early Extinction CS+ Fix Count (High IU Group)    1.09     1.81 0.40
## Early Extinction CS- Fix Count (High IU Group)    0.36    -0.42 0.39
## Late Extinction CS+ Fix Count (High IU Group)     1.00     1.61 0.44
## Late Extinction CS- Fix Count (High IU Group)     1.03     1.59 0.40
## Acquisition CS+ Fix Count (Low IU Group)          0.75    -0.01 0.40
## Acquisition CS- Fix Count (Low IU Group)          0.79    -0.13 0.39
## Early Extinction CS+ Fix Count (Low IU Group)     1.30     1.65 0.48
## Early Extinction CS- Fix Count (Low IU Group)     1.32     2.59 0.43
## Late Extinction CS+ Fix Count (Low IU Group)      0.47    -0.15 0.37
## Late Extinction CS- Fix Count (Low IU Group)      0.59    -0.05 0.41
# write to csv
write.csv(descriptives_fix_count_table, file = "tables/descriptives/descriptives_fix_count_table.csv",
          row.names = TRUE)

Fixation Duration

# for all participants
descriptives_all_fix_duration <- 
  describe(df[, c("acq_csp_fix_duration","acq_csm_fix_duration", 
                  "e_ext_csp_fix_duration", "e_ext_csm_fix_duration",
                  "l_ext_csp_fix_duration", "l_ext_csm_fix_duration")], 
           na.rm = TRUE)
  
# for high IU group
descriptives_high_iu_fix_duration <- 
  describe(df[df$iu_group == "1", c("acq_csp_fix_duration","acq_csm_fix_duration", 
                                    "e_ext_csp_fix_duration", "e_ext_csm_fix_duration", 
                                    "l_ext_csp_fix_duration", "l_ext_csm_fix_duration")], 
           na.rm = TRUE)
  
# for low IU group
descriptives_low_iu_fix_duration <- 
  describe(df[df$iu_group == "-1", c("acq_csp_fix_duration","acq_csm_fix_duration", 
                                     "e_ext_csp_fix_duration", "e_ext_csm_fix_duration",
                                     "l_ext_csp_fix_duration", "l_ext_csm_fix_duration")], 
           na.rm = TRUE)

# combine all in a table
descriptives_fix_duration_table <- round(rbind(descriptives_all_fix_duration, 
                                            descriptives_high_iu_fix_duration,
                                            descriptives_low_iu_fix_duration), 2)

# rename rows for easier interpretation
rownames(descriptives_fix_duration_table) <- c("Acquisition CS+ Fix Duration (All Participants)", 
                                            "Acquisition CS- Fix Duration (All Participants)", 
                                            "Early Extinction CS+ Fix Duration (All Participants)", 
                                            "Early Extinction CS- Fix Duration (All Participants)",
                                            "Late Extinction CS+ Fix Duration (All Participants)", 
                                            "Late Extinction CS- Fix Duration (All Participants)",
                                            "Acquisition CS+ Fix Duration (High IU Group)", 
                                            "Acquisition CS- Fix Duration (High IU Group)", 
                                            "Early Extinction CS+ Fix Duration (High IU Group)", 
                                            "Early Extinction CS- Fix Duration (High IU Group)",
                                            "Late Extinction CS+ Fix Duration (High IU Group)", 
                                            "Late Extinction CS- Fix Duration (High IU Group)",
                                            "Acquisition CS+ Fix Duration (Low IU Group)", 
                                            "Acquisition CS- Fix Duration (Low IU Group)", 
                                            "Early Extinction CS+ Fix Duration (Low IU Group)", 
                                            "Early Extinction CS- Fix Duration (Low IU Group)",
                                            "Late Extinction CS+ Fix Duration (Low IU Group)", 
                                            "Late Extinction CS- Fix Duration (Low IU Group)")

descriptives_fix_duration_table
##                                                      vars   n    mean      sd
## Acquisition CS+ Fix Duration (All Participants)         1 139 1309.36 1173.03
## Acquisition CS- Fix Duration (All Participants)         2 139 1200.18 1048.80
## Early Extinction CS+ Fix Duration (All Participants)    3 139 1104.04  930.02
## Early Extinction CS- Fix Duration (All Participants)    4 139 1203.66 1288.87
## Late Extinction CS+ Fix Duration (All Participants)     5 139 1066.13 1094.12
## Late Extinction CS- Fix Duration (All Participants)     6 139 1068.60 1204.27
## Acquisition CS+ Fix Duration (High IU Group)            1  68 1153.24 1126.41
## Acquisition CS- Fix Duration (High IU Group)            2  68 1003.87  938.91
## Early Extinction CS+ Fix Duration (High IU Group)       3  68  869.88  621.12
## Early Extinction CS- Fix Duration (High IU Group)       4  68  833.27  912.11
## Late Extinction CS+ Fix Duration (High IU Group)        5  68  799.03  732.10
## Late Extinction CS- Fix Duration (High IU Group)        6  68  719.91  687.99
## Acquisition CS+ Fix Duration (Low IU Group)             1  71 1458.89 1204.96
## Acquisition CS- Fix Duration (Low IU Group)             2  71 1388.19 1118.70
## Early Extinction CS+ Fix Duration (Low IU Group)        3  71 1328.31 1109.79
## Early Extinction CS- Fix Duration (Low IU Group)        4  71 1558.41 1489.19
## Late Extinction CS+ Fix Duration (Low IU Group)         5  71 1321.94 1308.18
## Late Extinction CS- Fix Duration (Low IU Group)         6  71 1402.56 1474.73
##                                                       median trimmed    mad
## Acquisition CS+ Fix Duration (All Participants)       789.44 1121.85 639.02
## Acquisition CS- Fix Duration (All Participants)       778.02 1017.24 606.82
## Early Extinction CS+ Fix Duration (All Participants)  786.25  958.21 611.53
## Early Extinction CS- Fix Duration (All Participants)  674.06  937.16 519.93
## Late Extinction CS+ Fix Duration (All Participants)   657.98  830.55 472.15
## Late Extinction CS- Fix Duration (All Participants)   578.00  786.22 397.29
## Acquisition CS+ Fix Duration (High IU Group)          666.40  943.03 510.01
## Acquisition CS- Fix Duration (High IU Group)          649.90  830.74 412.10
## Early Extinction CS+ Fix Duration (High IU Group)     716.65  780.80 516.83
## Early Extinction CS- Fix Duration (High IU Group)     533.31  641.57 308.46
## Late Extinction CS+ Fix Duration (High IU Group)      541.37  664.26 334.40
## Late Extinction CS- Fix Duration (High IU Group)      510.20  585.21 306.25
## Acquisition CS+ Fix Duration (Low IU Group)          1002.75 1309.89 888.32
## Acquisition CS- Fix Duration (Low IU Group)          1081.07 1216.87 977.62
## Early Extinction CS+ Fix Duration (Low IU Group)      931.86 1181.94 830.29
## Early Extinction CS- Fix Duration (Low IU Group)     1017.70 1282.48 964.06
## Late Extinction CS+ Fix Duration (Low IU Group)       781.46 1064.33 638.13
## Late Extinction CS- Fix Duration (Low IU Group)       845.97 1102.87 689.98
##                                                         min     max   range
## Acquisition CS+ Fix Duration (All Participants)       87.39 5083.82 4996.43
## Acquisition CS- Fix Duration (All Participants)       88.43 5446.78 5358.35
## Early Extinction CS+ Fix Duration (All Participants)  79.01 5346.50 5267.49
## Early Extinction CS- Fix Duration (All Participants)  65.23 6015.75 5950.52
## Late Extinction CS+ Fix Duration (All Participants)  121.30 5923.00 5801.70
## Late Extinction CS- Fix Duration (All Participants)  109.30 6086.56 5977.26
## Acquisition CS+ Fix Duration (High IU Group)          87.39 5083.82 4996.43
## Acquisition CS- Fix Duration (High IU Group)          88.43 5446.78 5358.35
## Early Extinction CS+ Fix Duration (High IU Group)    110.04 3044.00 2933.96
## Early Extinction CS- Fix Duration (High IU Group)     65.23 6015.75 5950.52
## Late Extinction CS+ Fix Duration (High IU Group)     121.84 4252.33 4130.49
## Late Extinction CS- Fix Duration (High IU Group)     109.30 4299.36 4190.06
## Acquisition CS+ Fix Duration (Low IU Group)          129.50 4985.33 4855.84
## Acquisition CS- Fix Duration (Low IU Group)          180.65 5219.17 5038.51
## Early Extinction CS+ Fix Duration (Low IU Group)      79.01 5346.50 5267.49
## Early Extinction CS- Fix Duration (Low IU Group)     119.97 5954.83 5834.86
## Late Extinction CS+ Fix Duration (Low IU Group)      121.30 5923.00 5801.70
## Late Extinction CS- Fix Duration (Low IU Group)      203.15 6086.56 5883.41
##                                                      skew kurtosis     se
## Acquisition CS+ Fix Duration (All Participants)      1.41     1.29  99.50
## Acquisition CS- Fix Duration (All Participants)      1.65     2.73  88.96
## Early Extinction CS+ Fix Duration (All Participants) 1.58     2.76  78.88
## Early Extinction CS- Fix Duration (All Participants) 2.05     3.83 109.32
## Late Extinction CS+ Fix Duration (All Participants)  2.17     4.52  92.80
## Late Extinction CS- Fix Duration (All Participants)  2.31     4.85 102.14
## Acquisition CS+ Fix Duration (High IU Group)         1.94     3.43 136.60
## Acquisition CS- Fix Duration (High IU Group)         2.30     6.41 113.86
## Early Extinction CS+ Fix Duration (High IU Group)    1.50     2.47  75.32
## Early Extinction CS- Fix Duration (High IU Group)    3.41    14.24 110.61
## Late Extinction CS+ Fix Duration (High IU Group)     2.66     8.49  88.78
## Late Extinction CS- Fix Duration (High IU Group)     3.05    10.94  83.43
## Acquisition CS+ Fix Duration (Low IU Group)          0.98    -0.04 143.00
## Acquisition CS- Fix Duration (Low IU Group)          1.21     1.01 132.77
## Early Extinction CS+ Fix Duration (Low IU Group)     1.16     1.00 131.71
## Early Extinction CS- Fix Duration (Low IU Group)     1.41     1.11 176.73
## Late Extinction CS+ Fix Duration (Low IU Group)      1.64     1.90 155.25
## Late Extinction CS- Fix Duration (Low IU Group)      1.64     1.62 175.02
# write to csv
write.csv(descriptives_fix_duration_table, file = "tables/descriptives/descriptives_fix_duration_table.csv", 
          row.names = TRUE)

Saccade Amplitude

# for all participants
descriptives_all_sacc_amplitude <- 
  describe(df[, c("acq_csp_sacc_amplitude","acq_csm_sacc_amplitude", 
                  "e_ext_csp_sacc_amplitude", "e_ext_csm_sacc_amplitude",
                  "l_ext_csp_sacc_amplitude", "l_ext_csm_sacc_amplitude")], 
           na.rm = TRUE)
  
# for high IU group
descriptives_high_iu_sacc_amplitude <- 
  describe(df[df$iu_group == "1", c("acq_csp_sacc_amplitude","acq_csm_sacc_amplitude",
                                    "e_ext_csp_sacc_amplitude", "e_ext_csm_sacc_amplitude",
                                    "l_ext_csp_sacc_amplitude", "l_ext_csm_sacc_amplitude")], 
           na.rm = TRUE)
  
# for low IU group
descriptives_low_iu_sacc_amplitude <- 
  describe(df[df$iu_group == "-1", c("acq_csp_sacc_amplitude","acq_csm_sacc_amplitude", 
                                     "e_ext_csp_sacc_amplitude", "e_ext_csm_sacc_amplitude",
                                     "l_ext_csp_sacc_amplitude", "l_ext_csm_sacc_amplitude")], 
           na.rm = TRUE)

# combine all into one table
descriptives_sacc_amplitude_table <- round(rbind(descriptives_all_sacc_amplitude, 
                                            descriptives_high_iu_sacc_amplitude,
                                            descriptives_low_iu_sacc_amplitude), 2)

# rename rows for easier interpretation
rownames(descriptives_sacc_amplitude_table) <- c("Acquisition CS+ Sacc Amplitude (All Participants)", 
                                            "Acquisition CS- Sacc Amplitude (All Participants)", 
                                            "Early Extinction CS+ Sacc Amplitude (All Participants)", 
                                            "Early Extinction CS- Sacc Amplitude (All Participants)",
                                            "Late Extinction CS+ Sacc Amplitude (All Participants)", 
                                            "Late Extinction CS- Sacc Amplitude (All Participants)",
                                            "Acquisition CS+ Sacc Amplitude (High IU Group)", 
                                            "Acquisition CS- Sacc Amplitude (High IU Group)", 
                                            "Early Extinction CS+ Sacc Amplitude (High IU Group)", 
                                            "Early Extinction CS- Sacc Amplitude (High IU Group)",
                                            "Late Extinction CS+ Sacc Amplitude (High IU Group)", 
                                            "Late Extinction CS- Sacc Amplitude (High IU Group)",
                                            "Acquisition CS+ Sacc Amplitude (Low IU Group)", 
                                            "Acquisition CS- Sacc Amplitude (Low IU Group)", 
                                            "Early Extinction CS+ Sacc Amplitude (Low IU Group)", 
                                            "Early Extinction CS- Sacc Amplitude (Low IU Group)",
                                            "Late Extinction CS+ Sacc Amplitude (Low IU Group)", 
                                            "Late Extinction CS- Sacc Amplitude (Low IU Group)")

descriptives_sacc_amplitude_table
##                                                        vars   n mean   sd
## Acquisition CS+ Sacc Amplitude (All Participants)         1 139 2.88 1.51
## Acquisition CS- Sacc Amplitude (All Participants)         2 137 2.98 1.57
## Early Extinction CS+ Sacc Amplitude (All Participants)    3 138 3.07 1.81
## Early Extinction CS- Sacc Amplitude (All Participants)    4 138 3.13 1.73
## Late Extinction CS+ Sacc Amplitude (All Participants)     5 138 3.00 1.92
## Late Extinction CS- Sacc Amplitude (All Participants)     6 139 3.10 1.97
## Acquisition CS+ Sacc Amplitude (High IU Group)            1  68 3.10 1.71
## Acquisition CS- Sacc Amplitude (High IU Group)            2  67 3.16 1.70
## Early Extinction CS+ Sacc Amplitude (High IU Group)       3  68 3.21 1.80
## Early Extinction CS- Sacc Amplitude (High IU Group)       4  67 3.46 1.88
## Late Extinction CS+ Sacc Amplitude (High IU Group)        5  68 3.21 1.78
## Late Extinction CS- Sacc Amplitude (High IU Group)        6  68 3.37 1.85
## Acquisition CS+ Sacc Amplitude (Low IU Group)             1  71 2.66 1.27
## Acquisition CS- Sacc Amplitude (Low IU Group)             2  70 2.80 1.43
## Early Extinction CS+ Sacc Amplitude (Low IU Group)        3  70 2.95 1.83
## Early Extinction CS- Sacc Amplitude (Low IU Group)        4  71 2.81 1.53
## Late Extinction CS+ Sacc Amplitude (Low IU Group)         5  70 2.79 2.03
## Late Extinction CS- Sacc Amplitude (Low IU Group)         6  71 2.84 2.06
##                                                        median trimmed  mad  min
## Acquisition CS+ Sacc Amplitude (All Participants)        2.64    2.71 1.35 0.43
## Acquisition CS- Sacc Amplitude (All Participants)        2.65    2.81 1.25 0.54
## Early Extinction CS+ Sacc Amplitude (All Participants)   2.78    2.85 1.65 0.58
## Early Extinction CS- Sacc Amplitude (All Participants)   2.66    2.94 1.35 0.42
## Late Extinction CS+ Sacc Amplitude (All Participants)    2.69    2.78 1.61 0.38
## Late Extinction CS- Sacc Amplitude (All Participants)    2.90    2.86 1.89 0.42
## Acquisition CS+ Sacc Amplitude (High IU Group)           2.99    2.92 1.74 0.43
## Acquisition CS- Sacc Amplitude (High IU Group)           2.86    2.96 1.49 0.54
## Early Extinction CS+ Sacc Amplitude (High IU Group)      3.08    3.02 1.75 0.64
## Early Extinction CS- Sacc Amplitude (High IU Group)      3.18    3.28 1.74 0.69
## Late Extinction CS+ Sacc Amplitude (High IU Group)       2.90    3.03 1.78 0.38
## Late Extinction CS- Sacc Amplitude (High IU Group)       3.13    3.22 1.90 0.61
## Acquisition CS+ Sacc Amplitude (Low IU Group)            2.52    2.56 1.17 0.59
## Acquisition CS- Sacc Amplitude (Low IU Group)            2.60    2.66 1.23 0.59
## Early Extinction CS+ Sacc Amplitude (Low IU Group)       2.48    2.70 1.52 0.58
## Early Extinction CS- Sacc Amplitude (Low IU Group)       2.34    2.65 1.07 0.42
## Late Extinction CS+ Sacc Amplitude (Low IU Group)        2.25    2.52 1.69 0.43
## Late Extinction CS- Sacc Amplitude (Low IU Group)        2.63    2.53 1.93 0.42
##                                                          max range skew
## Acquisition CS+ Sacc Amplitude (All Participants)       8.15  7.72 1.01
## Acquisition CS- Sacc Amplitude (All Participants)       8.57  8.04 1.12
## Early Extinction CS+ Sacc Amplitude (All Participants)  9.18  8.59 1.14
## Early Extinction CS- Sacc Amplitude (All Participants) 11.42 11.00 1.44
## Late Extinction CS+ Sacc Amplitude (All Participants)  13.11 12.73 1.72
## Late Extinction CS- Sacc Amplitude (All Participants)  10.95 10.53 1.25
## Acquisition CS+ Sacc Amplitude (High IU Group)          8.15  7.72 0.96
## Acquisition CS- Sacc Amplitude (High IU Group)          8.57  8.04 1.14
## Early Extinction CS+ Sacc Amplitude (High IU Group)     8.65  8.01 0.89
## Early Extinction CS- Sacc Amplitude (High IU Group)    11.42 10.73 1.38
## Late Extinction CS+ Sacc Amplitude (High IU Group)      9.62  9.24 1.11
## Late Extinction CS- Sacc Amplitude (High IU Group)      9.74  9.13 0.93
## Acquisition CS+ Sacc Amplitude (Low IU Group)           6.35  5.76 0.68
## Acquisition CS- Sacc Amplitude (Low IU Group)           7.37  6.78 0.93
## Early Extinction CS+ Sacc Amplitude (Low IU Group)      9.18  8.59 1.37
## Early Extinction CS- Sacc Amplitude (Low IU Group)      9.11  8.69 1.33
## Late Extinction CS+ Sacc Amplitude (Low IU Group)      13.11 12.68 2.20
## Late Extinction CS- Sacc Amplitude (Low IU Group)      10.95 10.53 1.57
##                                                        kurtosis   se
## Acquisition CS+ Sacc Amplitude (All Participants)          0.86 0.13
## Acquisition CS- Sacc Amplitude (All Participants)          1.35 0.13
## Early Extinction CS+ Sacc Amplitude (All Participants)     1.13 0.15
## Early Extinction CS- Sacc Amplitude (All Participants)     3.30 0.15
## Late Extinction CS+ Sacc Amplitude (All Participants)      5.31 0.16
## Late Extinction CS- Sacc Amplitude (All Participants)      1.95 0.17
## Acquisition CS+ Sacc Amplitude (High IU Group)             0.40 0.21
## Acquisition CS- Sacc Amplitude (High IU Group)             1.13 0.21
## Early Extinction CS+ Sacc Amplitude (High IU Group)        0.32 0.22
## Early Extinction CS- Sacc Amplitude (High IU Group)        3.07 0.23
## Late Extinction CS+ Sacc Amplitude (High IU Group)         1.55 0.22
## Late Extinction CS- Sacc Amplitude (High IU Group)         0.76 0.22
## Acquisition CS+ Sacc Amplitude (Low IU Group)             -0.16 0.15
## Acquisition CS- Sacc Amplitude (Low IU Group)              0.81 0.17
## Early Extinction CS+ Sacc Amplitude (Low IU Group)         1.92 0.22
## Early Extinction CS- Sacc Amplitude (Low IU Group)         2.40 0.18
## Late Extinction CS+ Sacc Amplitude (Low IU Group)          7.87 0.24
## Late Extinction CS- Sacc Amplitude (Low IU Group)          3.05 0.25
# write to csv
write.csv(descriptives_sacc_amplitude_table, file = "tables/descriptives/descriptives_sacc_amplitude_table.csv", 
          row.names = TRUE)

Data Transformation

Log-Transformation of Fixation Duration

# as fixation duration had high skew (>3) in high IU group for early and late 
# extinction CS-, fixation duration will be log-transformed for each condition

# for acquisition CS+
df$acq_csp_fix_duration_log <- log(df$acq_csp_fix_duration)

# for acquisition CS-
df$acq_csm_fix_duration_log <- log(df$acq_csm_fix_duration)

# for early extinction CS+
df$e_ext_csp_fix_duration_log <- log(df$e_ext_csp_fix_duration)

# for early extinction CS-
df$e_ext_csm_fix_duration_log <- log(df$e_ext_csm_fix_duration)

# for late extinction CS+
df$l_ext_csp_fix_duration_log <- log(df$l_ext_csp_fix_duration)

# for late extinction CS-
df$l_ext_csm_fix_duration_log <- log(df$l_ext_csm_fix_duration)

Check Descriptives of Fixation Duration Following Log-Transformation

# re-compute descriptives for fixation duration following log transformation

# for all participants
descriptives_all_fix_duration_log <- 
  describe(df[, c("acq_csp_fix_duration_log","acq_csm_fix_duration_log", 
                  "e_ext_csp_fix_duration_log", "e_ext_csm_fix_duration_log",
                  "l_ext_csp_fix_duration_log", "l_ext_csm_fix_duration_log")], 
           na.rm = TRUE)
  
# for high IU group
descriptives_high_iu_fix_duration_log <- 
  describe(df[df$iu_group == "1", c("acq_csp_fix_duration_log","acq_csm_fix_duration_log", 
                                    "e_ext_csp_fix_duration_log", "e_ext_csm_fix_duration_log", 
                                    "l_ext_csp_fix_duration_log", "l_ext_csm_fix_duration_log")], 
           na.rm = TRUE)
  
# for low IU group
descriptives_low_iu_fix_duration_log <- 
  describe(df[df$iu_group == "-1", c("acq_csp_fix_duration_log","acq_csm_fix_duration_log", 
                                     "e_ext_csp_fix_duration_log", "e_ext_csm_fix_duration_log",
                                     "l_ext_csp_fix_duration_log", "l_ext_csm_fix_duration_log")], 
           na.rm = TRUE)

# combine all to table
descriptives_fix_duration_table_log <- round(rbind(descriptives_all_fix_duration_log, 
                                            descriptives_high_iu_fix_duration_log,
                                            descriptives_low_iu_fix_duration_log), 2)

# rename rows for easier interpretation
rownames(descriptives_fix_duration_table_log) <- c("Acquisition CS+ Fix Duration (All Participants)", 
                                            "Acquisition CS- Fix Duration (All Participants)", 
                                            "Early Extinction CS+ Fix Duration (All Participants)", 
                                            "Early Extinction CS- Fix Duration (All Participants)",
                                            "Late Extinction CS+ Fix Duration (All Participants)", 
                                            "Late Extinction CS- Fix Duration (All Participants)",
                                            "Acquisition CS+ Fix Duration (High IU Group)", 
                                            "Acquisition CS- Fix Duration (High IU Group)", 
                                            "Early Extinction CS+ Fix Duration (High IU Group)", 
                                            "Early Extinction CS- Fix Duration (High IU Group)",
                                            "Late Extinction CS+ Fix Duration (High IU Group)", 
                                            "Late Extinction CS- Fix Duration (High IU Group)",
                                            "Acquisition CS+ Fix Duration (Low IU Group)", 
                                            "Acquisition CS- Fix Duration (Low IU Group)", 
                                            "Early Extinction CS+ Fix Duration (Low IU Group)", 
                                            "Early Extinction CS- Fix Duration (Low IU Group)",
                                            "Late Extinction CS+ Fix Duration (Low IU Group)", 
                                            "Late Extinction CS- Fix Duration (Low IU Group)")

descriptives_fix_duration_table_log
##                                                      vars   n mean   sd median
## Acquisition CS+ Fix Duration (All Participants)         1 139 6.80 0.89   6.67
## Acquisition CS- Fix Duration (All Participants)         2 139 6.76 0.83   6.66
## Early Extinction CS+ Fix Duration (All Participants)    3 139 6.68 0.84   6.67
## Early Extinction CS- Fix Duration (All Participants)    4 139 6.67 0.90   6.51
## Late Extinction CS+ Fix Duration (All Participants)     5 139 6.60 0.83   6.49
## Late Extinction CS- Fix Duration (All Participants)     6 139 6.56 0.85   6.36
## Acquisition CS+ Fix Duration (High IU Group)            1  68 6.68 0.87   6.50
## Acquisition CS- Fix Duration (High IU Group)            2  68 6.60 0.78   6.48
## Early Extinction CS+ Fix Duration (High IU Group)       3  68 6.54 0.70   6.57
## Early Extinction CS- Fix Duration (High IU Group)       4  68 6.40 0.75   6.28
## Late Extinction CS+ Fix Duration (High IU Group)        5  68 6.41 0.72   6.29
## Late Extinction CS- Fix Duration (High IU Group)        6  68 6.31 0.71   6.23
## Acquisition CS+ Fix Duration (Low IU Group)             1  71 6.92 0.89   6.91
## Acquisition CS- Fix Duration (Low IU Group)             2  71 6.91 0.85   6.99
## Early Extinction CS+ Fix Duration (Low IU Group)        3  71 6.81 0.94   6.84
## Early Extinction CS- Fix Duration (Low IU Group)        4  71 6.92 0.96   6.93
## Late Extinction CS+ Fix Duration (Low IU Group)         5  71 6.79 0.89   6.66
## Late Extinction CS- Fix Duration (Low IU Group)         6  71 6.81 0.91   6.74
##                                                      trimmed  mad  min  max
## Acquisition CS+ Fix Duration (All Participants)         6.80 0.94 4.47 8.53
## Acquisition CS- Fix Duration (All Participants)         6.74 0.97 4.48 8.60
## Early Extinction CS+ Fix Duration (All Participants)    6.69 0.93 4.37 8.58
## Early Extinction CS- Fix Duration (All Participants)    6.63 0.84 4.18 8.70
## Late Extinction CS+ Fix Duration (All Participants)     6.55 0.86 4.80 8.69
## Late Extinction CS- Fix Duration (All Participants)     6.50 0.71 4.69 8.71
## Acquisition CS+ Fix Duration (High IU Group)            6.67 0.79 4.47 8.53
## Acquisition CS- Fix Duration (High IU Group)            6.57 0.74 4.48 8.60
## Early Extinction CS+ Fix Duration (High IU Group)       6.54 0.77 4.70 8.02
## Early Extinction CS- Fix Duration (High IU Group)       6.35 0.67 4.18 8.70
## Late Extinction CS+ Fix Duration (High IU Group)        6.38 0.75 4.80 8.36
## Late Extinction CS- Fix Duration (High IU Group)        6.29 0.61 4.69 8.37
## Acquisition CS+ Fix Duration (Low IU Group)             6.94 1.15 4.86 8.51
## Acquisition CS- Fix Duration (Low IU Group)             6.91 1.01 5.20 8.56
## Early Extinction CS+ Fix Duration (Low IU Group)        6.86 1.15 4.37 8.58
## Early Extinction CS- Fix Duration (Low IU Group)        6.91 1.14 4.79 8.69
## Late Extinction CS+ Fix Duration (Low IU Group)         6.75 0.91 4.80 8.69
## Late Extinction CS- Fix Duration (Low IU Group)         6.75 0.96 5.31 8.71
##                                                      range  skew kurtosis   se
## Acquisition CS+ Fix Duration (All Participants)       4.06  0.02    -0.64 0.08
## Acquisition CS- Fix Duration (All Participants)       4.12  0.11    -0.65 0.07
## Early Extinction CS+ Fix Duration (All Participants)  4.21 -0.11    -0.40 0.07
## Early Extinction CS- Fix Duration (All Participants)  4.52  0.34    -0.31 0.08
## Late Extinction CS+ Fix Duration (All Participants)   3.89  0.45    -0.27 0.07
## Late Extinction CS- Fix Duration (All Participants)   4.02  0.59    -0.02 0.07
## Acquisition CS+ Fix Duration (High IU Group)          4.06  0.09    -0.21 0.11
## Acquisition CS- Fix Duration (High IU Group)          4.12  0.25    -0.01 0.09
## Early Extinction CS+ Fix Duration (High IU Group)     3.32 -0.08    -0.48 0.08
## Early Extinction CS- Fix Duration (High IU Group)     4.52  0.48     1.07 0.09
## Late Extinction CS+ Fix Duration (High IU Group)      3.55  0.41    -0.03 0.09
## Late Extinction CS- Fix Duration (High IU Group)      3.67  0.35     0.66 0.09
## Acquisition CS+ Fix Duration (Low IU Group)           3.65 -0.07    -1.03 0.11
## Acquisition CS- Fix Duration (Low IU Group)           3.36 -0.07    -1.05 0.10
## Early Extinction CS+ Fix Duration (Low IU Group)      4.21 -0.32    -0.53 0.11
## Early Extinction CS- Fix Duration (Low IU Group)      3.90  0.03    -0.87 0.11
## Late Extinction CS+ Fix Duration (Low IU Group)       3.89  0.29    -0.67 0.11
## Late Extinction CS- Fix Duration (Low IU Group)       3.40  0.47    -0.83 0.11
# write to csv
write.csv(descriptives_fix_duration_table_log, file = "tables/descriptives/descriptives_fix_duration_table_log.csv", 
          row.names = TRUE)

### there are no longer any skew values of +/- 3. 

Check Histograms of Fixation Duration Following Log-Transformation

Acquisition CS+

########## pre-log-transformation
hist_acq_csp_fix_duration

########## post-log-transformation
hist_acq_csp_fix_duration_log <- df %>%
  ggplot(aes(acq_csp_fix_duration_log, fill = iu_group)) +
  geom_histogram(binwidth = .2, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 12, 2)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Acquisition CS+ Fixation Duration (Log-Transformed)") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_acq_csp_fix_duration_log

# save plot to file
ggsave(filename = "graphs/histograms/hist_acq_csp_fix_duration_log.png", 
       plot = hist_acq_csp_fix_duration_log, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Acquisition CS-

########## pre-log-transformation
hist_acq_csp_fix_duration_log

########## post-log-transformation
hist_acq_csm_fix_duration_log <- df %>%
  ggplot(aes(acq_csm_fix_duration_log, fill = iu_group)) +
  geom_histogram(binwidth = .2, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 12, 2)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Acquisition CS- Fixation Duration (Log-Transformed)") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_acq_csm_fix_duration_log

# save plot to file
ggsave(filename = "graphs/histograms/hist_acq_csm_fix_duration_log.png", 
       plot = hist_acq_csm_fix_duration_log, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

# combine histograms of acquisition fix duration pre and post log-transformation
hists_fix_duration_acq_log <- grid.arrange(hist_acq_csp_fix_duration, hist_acq_csm_fix_duration,
                                       hist_acq_csp_fix_duration_log, hist_acq_csm_fix_duration_log,
                                       ncol = 2)

# save plot to file
ggsave(filename = "graphs/histograms/hists_fix_duration_acq_log.png", 
       plot = hists_fix_duration_acq_log, 
       width = 30, 
       height = 20, 
       dpi = 300, 
       units = "cm")

Early Extinction CS+

########## pre-log-transformation
hist_e_ext_csp_fix_duration

########## post-log-transformation
hist_e_ext_csp_fix_duration_log <- df %>%
  ggplot(aes(e_ext_csp_fix_duration_log, fill = iu_group)) +
  geom_histogram(binwidth = .2, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 12, 2)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Early Extinction CS+ Fixation Duration (Log-Transformed)") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_e_ext_csp_fix_duration_log

# save plot to file
ggsave(filename = "graphs/histograms/hist_e_ext_csp_fix_duration_log.png", 
       plot = hist_e_ext_csp_fix_duration_log, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Early Extinction CS-

########## pre-log-transformation
hist_e_ext_csm_fix_duration

########## post-log-transformation
hist_e_ext_csm_fix_duration_log <- df %>%
  ggplot(aes(e_ext_csm_fix_duration_log, fill = iu_group)) +
  geom_histogram(binwidth = .2, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 12, 2)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Early Extinction CS- Fixation Duration (Log-Transformed)") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_e_ext_csm_fix_duration_log

# save plot to file
ggsave(filename = "graphs/histograms/hist_e_ext_csm_fix_duration_log.png", 
       plot = hist_e_ext_csm_fix_duration_log, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

# combine histograms of early extinction fix duration pre and post log-transformation
hists_fix_duration_e_ext_log <- grid.arrange(hist_e_ext_csp_fix_duration, hist_e_ext_csm_fix_duration,
                                       hist_e_ext_csp_fix_duration_log, hist_e_ext_csm_fix_duration_log,
                                       ncol = 2)

# save plot to file
ggsave(filename = "graphs/histograms/hists_fix_duration_e_ext_log.png", 
       plot = hists_fix_duration_e_ext_log, 
       width = 30, 
       height = 20, 
       dpi = 300, 
       units = "cm")

Late Extinction CS+

########## pre-log-transformation
hist_l_ext_csp_fix_duration

########## post-log-transformation
hist_l_ext_csp_fix_duration_log <- df %>%
  ggplot(aes(l_ext_csp_fix_duration_log, fill = iu_group)) +
  geom_histogram(binwidth = .2, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 12, 2)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Late Extinction CS+ Fixation Duration (Log-Transformed)") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_l_ext_csp_fix_duration_log

# save plot to file
ggsave(filename = "graphs/histograms/hist_l_ext_csp_fix_duration_log.png", 
       plot = hist_l_ext_csp_fix_duration_log, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Late Extinction CS-

########## pre-log-transformation
hist_l_ext_csm_fix_duration

########## post-log-transformation
hist_l_ext_csm_fix_duration_log <- df %>%
  ggplot(aes(l_ext_csm_fix_duration_log, fill = iu_group)) +
  geom_histogram(binwidth = .2, colour = "white", alpha = .5, position = "identity") +
  theme_classic() +
  theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  scale_x_continuous(breaks = seq(0, 12, 2)) +
  labs(x = "Fixation Duration (ms)", y = "Frequency") +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  ggtitle("Histogram Late Extinction CS- Fixation Duration (Log-Transformed)") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  guides(fill = guide_legend(reverse = TRUE)) +
  scale_fill_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
  labs(fill = "IU Group") 
hist_l_ext_csm_fix_duration_log

# save plot to file
ggsave(filename = "graphs/histograms/hist_l_ext_csm_fix_duration_log.png", 
       plot = hist_l_ext_csm_fix_duration_log, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

# combine histograms of late extinction fix duration pre and post log-transformation
hists_fix_duration_l_ext_log <- grid.arrange(hist_l_ext_csp_fix_duration, hist_l_ext_csm_fix_duration,
                                       hist_l_ext_csp_fix_duration_log, hist_l_ext_csm_fix_duration_log,
                                       ncol = 2)

# save plot to file
ggsave(filename = "graphs/histograms/hists_fix_duration_l_ext_log.png", 
       plot = hists_fix_duration_l_ext_log, 
       width = 30, 
       height = 20, 
       dpi = 300, 
       units = "cm")

ANOVAs

ANOVA Acquisition Fixation Count

# transform wide format data into long format for mixed ANOVA 
df_long_acq_fix_count <- melt(df, id = c("id", "iu_group"), 
                                 measure.vars = c("acq_csp_fix_count", 
                                                  "acq_csm_fix_count"))

# rename columns for easier interpretation
colnames(df_long_acq_fix_count) = c("id", "iu_group", "condition", "fix_count")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_acq_fix_count$stimulus <- 
  factor(ifelse(df_long_acq_fix_count$condition == "acq_csp_fix_count", 1, -1))

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) mixed ANOVA, 
# and obtain effect size (partial eta squared)
acq_fix_count_anova <- 
  anova_test(df_long_acq_fix_count, fix_count ~ iu_group * stimulus + Error(id/stimulus),
                        effect.size = "pes") 

# obtain the mixed ANOVA results
get_anova_table(acq_fix_count_anova)
## ANOVA Table (type III tests)
## 
##              Effect DFn DFd      F        p p<.05   pes
## 1          iu_group   1 137  4.806 0.030000     * 0.034
## 2          stimulus   1 137 11.441 0.000937     * 0.077
## 3 iu_group:stimulus   1 137  0.258 0.613000       0.002
# results:
# IU: F(1,137) = 4.81, p = .030*, eta2(partial) = .034
# Stimulus: F(1,137) = 11.44, p < .001***, eta2(partial) = .077
# IU * Stimulus: F(1, 137) = 0.26, p = .613, eta2(partial) = .002

# therefore, there is a significant effect of IU & Stimulus on fixation count in acquisition, 
# and no significant IU*Stimulus interaction

# write to csv
write.csv((get_anova_table(acq_fix_count_anova)), 
          file = "tables/anovas/acq_fix_count_anova.csv")

ANOVA Acquisition Fixation Duration (Log Transformed)

# transform wide format data into long format for mixed ANOVA 
df_long_acq_fix_duration_log <- melt(df, id = c("id", "iu_group"), 
                                 measure.vars = c("acq_csp_fix_duration_log", 
                                                  "acq_csm_fix_duration_log"))

# the first input for melt command is the df with wide format. Second input
# is id =, which is where we list ppts with specific variables within wide format
# df. Here we have ppts ID no's as participant specific variable, and IU 
# group they are assigned to is also specific for each participant. Therefore
# only need to list the two variables after id =. 

# rename columns for easier interpretation
colnames(df_long_acq_fix_duration_log) = c("id", "iu_group", "condition", "fix_duration_log")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_acq_fix_duration_log$stimulus <- 
  factor(ifelse(df_long_acq_fix_duration_log$condition == "acq_csp_fix_duration_log", 1, -1))

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) mixed ANOVA, 
# and obtain effect size (partial eta squared)
acq_fix_duration_anova_log <- 
  anova_test(df_long_acq_fix_duration_log, fix_duration_log ~ iu_group * stimulus + Error(id/stimulus),
                        effect.size = "pes") 
# the error(id/stimulus) variable is unique to repeated-measures ANOVA, and means
# that the variable 'stimulus' is manipulated within 'id'

# obtain the mixed ANOVA results
get_anova_table(acq_fix_duration_anova_log)
## ANOVA Table (type III tests)
## 
##              Effect DFn DFd     F     p p<.05   pes
## 1          iu_group   1 137 3.907 0.050       0.028
## 2          stimulus   1 137 2.921 0.090       0.021
## 3 iu_group:stimulus   1 137 1.271 0.261       0.009
# results:
# IU: F(1,137) = 3.91, p = .050*, eta2(partial) = .028
# Stimulus: F(1,137) = 2.92, p = .090, eta2(partial) = .021
# IU * Stimulus: F(1, 137) = 1.27, p = .261, eta2(partial) = .009

# therefore, there is a sig effect of IU, and no
# sig effect of stimulus or IU-stimulus interaction 

# write to csv
write.csv((get_anova_table(acq_fix_duration_anova_log)), 
          file = "tables/anovas/acq_fix_duration_anova_log.csv")

ANOVA Acquisition Saccade Amplitude

# transform wide format data into long format for mixed ANOVA 
df_long_acq_sacc_amplitude <- melt(df, id = c("id", "iu_group"), 
                                 measure.vars = c("acq_csp_sacc_amplitude", 
                                                  "acq_csm_sacc_amplitude"))

# rename columns for easier interpretation
colnames(df_long_acq_sacc_amplitude) = c("id", "iu_group", "condition", "sacc_amplitude")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_acq_sacc_amplitude$stimulus <- 
  factor(ifelse(df_long_acq_sacc_amplitude$condition == "acq_csp_sacc_amplitude", 1, -1))

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) mixed ANOVA, 
# and obtain effect size (partial eta squared)
acq_sacc_amplitude_anova <- 
  anova_test(df_long_acq_sacc_amplitude, sacc_amplitude ~ iu_group * stimulus + Error(id/stimulus),
             effect.size = "pes") 
## Warning: NA detected in rows: 234,259.
## Removing this rows before the analysis.
# obtain the mixed ANOVA results
get_anova_table(acq_sacc_amplitude_anova)
## ANOVA Table (type III tests)
## 
##              Effect DFn DFd     F     p p<.05   pes
## 1          iu_group   1 135 2.984 0.086       0.022
## 2          stimulus   1 135 0.950 0.332       0.007
## 3 iu_group:stimulus   1 135 0.379 0.539       0.003
# results:
# IU: F(1,135) = 2.98, p = .086, eta2(partial) = .022
# Stimulus: F(1,135) = 0.95, p = .332, eta2(partial) = .007
# IU * Stimulus: F(1, 135) = 0.38, p = .539, eta2(partial) = .003

# therefore, there are no significant effects on saccade amplitude in
# acquisition

# write to csv
write.csv((get_anova_table(acq_sacc_amplitude_anova)), 
          file = "tables/anovas/acq_sacc_amplitude_anova.csv")

ANOVA Extinction Fixation Count

# transform wide format data into long format for mixed ANOVA 
df_long_ext_fix_count <- melt(df, id = c("id", "iu_group"), 
                                 measure.vars = c("e_ext_csp_fix_count", 
                                                  "e_ext_csm_fix_count",
                                                  "l_ext_csp_fix_count", 
                                                  "l_ext_csm_fix_count"))

# rename columns for easier interpretation
colnames(df_long_ext_fix_count) = c("id", "iu_group", "condition", "fix_count")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_ext_fix_count$stimulus <- 
  factor(ifelse(df_long_ext_fix_count$condition == "e_ext_csp_fix_count" | 
                  df_long_ext_fix_count$condition == "l_ext_csp_fix_count", 1, -1))

# create column to code extinction as early (1) and late (-1)
df_long_ext_fix_count$time <- 
  factor(ifelse(df_long_ext_fix_count$condition == "e_ext_csp_fix_count" | 
                  df_long_ext_fix_count$condition == "e_ext_csm_fix_count", 1, -1))

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) x 2 (Time: Early, Late) mixed ANOVA, 
# and obtain effect size (partial eta squared)
ext_fix_count_anova <- 
  anova_test(df_long_ext_fix_count, 
             fix_count ~ iu_group * stimulus * time + Error(id/(stimulus*time)),
             effect.size = "pes") 

# obtain the mixed ANOVA results
get_anova_table(ext_fix_count_anova)
## ANOVA Table (type III tests)
## 
##                   Effect DFn DFd     F     p p<.05      pes
## 1               iu_group   1 137 7.672 0.006     * 0.053000
## 2               stimulus   1 137 4.155 0.043     * 0.029000
## 3                   time   1 137 5.733 0.018     * 0.040000
## 4      iu_group:stimulus   1 137 3.460 0.065       0.025000
## 5          iu_group:time   1 137 4.572 0.034     * 0.032000
## 6          stimulus:time   1 137 0.061 0.806       0.000443
## 7 iu_group:stimulus:time   1 137 0.600 0.440       0.004000
# results:
# IU: F(1,137) = 7.67, p = .006 ***, eta2(partial) = .053
# Stimulus: F(1,137) = 4.16, p = .043 *, eta2(partial) = .029
# Time: F(1,137) = 5.73, p = .018 *, eta2(partial) = .049
# IU * Stimulus: F(1, 137) = 3.46, p = .065, eta2(partial) = .025
# IU * Time: F(1,137) = 4.57, p = .034 *, eta2(partial) = . 032
# Stimulus * Time: F(1,137) = 0.06, p = .806, eta2(partial) < .001
# IU * Stimulus * Time: F(1,137) = 0.60, p = .440, eta2(partial) = .004

# therefore, there is a significant effect of IU, Stimulus and Time on fixation count in extinction, 
# as well as a significant interaction effect of IU * Time,
# but no other significant interactions. 

# write to csv
write.csv((get_anova_table(ext_fix_count_anova)), 
          file = "tables/anovas/ext_fix_count_anova.csv")

# as there was a significant IU*Time interaction, conduct simple 
# main effects analysis:
## obtain effect of IU at each level of time 
simple_effects_ext_fix_count_iu <- df_long_ext_fix_count %>%
  group_by(time) %>%
  anova_test(dv = fix_count, wid = id, between = iu_group, within = stimulus, effect.size = "pes") %>%
  get_anova_table() 

# get the output
simple_effects_ext_fix_count_iu
## # A tibble: 6 × 8
##   time  Effect              DFn   DFd      F        p `p<.05`   pes
## * <fct> <chr>             <dbl> <dbl>  <dbl>    <dbl> <chr>   <dbl>
## 1 -1    iu_group              1   137 11.4   0.000952 "*"     0.077
## 2 -1    stimulus              1   137  3.38  0.068    ""      0.024
## 3 -1    iu_group:stimulus     1   137  0.864 0.354    ""      0.006
## 4 1     iu_group              1   137  3.63  0.059    ""      0.026
## 5 1     stimulus              1   137  1.50  0.222    ""      0.011
## 6 1     iu_group:stimulus     1   137  3.04  0.084    ""      0.022
# results: 
# The effect of IU group at early extinction was not significant [F(1,137) = 3.63, p = .059, pes = .026]
# the effect of IU group at late extinction was significant [F(1,137) = 11.41, p < .001, pes = .077]

ANOVA Extinction Fixation Duration (Log Transformed)

# transform wide format data into long format for mixed ANOVA 
df_long_ext_fix_duration_log <- melt(df, id = c("id", "iu_group"), 
                                 measure.vars = c("e_ext_csp_fix_duration_log", 
                                                  "e_ext_csm_fix_duration_log",
                                                  "l_ext_csp_fix_duration_log", 
                                                  "l_ext_csm_fix_duration_log"))

# rename columns for easier interpretation
colnames(df_long_ext_fix_duration_log) = c("id", "iu_group", "condition", "fix_duration_log")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_ext_fix_duration_log$stimulus <- 
  factor(ifelse(df_long_ext_fix_duration_log$condition == "e_ext_csp_fix_duration_log" | 
                  df_long_ext_fix_duration_log$condition == "l_ext_csp_fix_duration_log", 1, -1))

# create column to code extinction as early (1) and late (-1)
df_long_ext_fix_duration_log$time <- 
  factor(ifelse(df_long_ext_fix_duration_log$condition == "e_ext_csp_fix_duration_log" | 
                  df_long_ext_fix_duration_log$condition == "e_ext_csm_fix_duration_log", 1, -1))

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) x 2 (Time: Early, Late) mixed ANOVA, 
# and obtain effect size (partial eta squared)
ext_fix_duration_anova_log <- 
  anova_test(df_long_ext_fix_duration_log, 
             fix_duration_log ~ iu_group * stimulus * time + Error(id/(stimulus*time)),
             effect.size = "pes") 

# obtain the mixed ANOVA results
get_anova_table(ext_fix_duration_anova_log)
## ANOVA Table (type III tests)
## 
##                   Effect DFn DFd      F     p p<.05   pes
## 1               iu_group   1 137 11.213 0.001     * 0.076
## 2               stimulus   1 137  0.510 0.477       0.004
## 3                   time   1 137  4.351 0.039     * 0.031
## 4      iu_group:stimulus   1 137  5.823 0.017     * 0.041
## 5          iu_group:time   1 137  0.241 0.624       0.002
## 6          stimulus:time   1 137  0.171 0.680       0.001
## 7 iu_group:stimulus:time   1 137  0.946 0.333       0.007
# results:
# IU: F(1,137) = 11.21, p < .001 *, eta2(partial) = .076
# Stimulus: F(1,137) = 0.51, p = .477, eta2(partial) = .004
# Time: F(1,137) = 4.35, p = .039*, eta2(partial) = .031
# IU * Stimulus: F(1, 137) = 5.82, p = .017*, eta2(partial) = .041
# IU * Time: F(1,137) = 0.24, p = .624, eta2(partial) = .002
# Stimulus * Time: F(1,137) = 0.17, p = 680, eta2(partial) = .001
# IU * Stimulus * Time: F(1,137) = 0.95, p = .333, eta2(partial) = .007

# therefore, there is a significant effect of IU, Time and IU-Stimulus
# interaction on fixation duration in extinction, 
# and no other significant effects or interactions. 

# write to csv
write.csv((get_anova_table(ext_fix_duration_anova_log)), 
          file = "tables/anovas/ext_fix_duration_anova_log.csv")

# as there was a significant IU*Stimulus interaction, conduct simple 
# main effects analysis:
## obtain effect of IU at each level of stimulus 
simple_effects_ext_fix_duration_log_iu <- df_long_ext_fix_duration_log %>%
  group_by(stimulus) %>%
  anova_test(dv = fix_duration_log, wid = id, between = iu_group, within = time, effect.size = "pes") %>%
  get_anova_table() 

# get the output
simple_effects_ext_fix_duration_log_iu
## # A tibble: 6 × 8
##   stimulus Effect          DFn   DFd      F        p `p<.05`      pes
## * <fct>    <chr>         <dbl> <dbl>  <dbl>    <dbl> <chr>      <dbl>
## 1 -1       iu_group          1   137 14.4   0.000218 "*"     0.095   
## 2 -1       time              1   137  4.34  0.039    "*"     0.031   
## 3 -1       iu_group:time     1   137  0.026 0.871    ""      0.000192
## 4 1        iu_group          1   137  6.70  0.011    "*"     0.047   
## 5 1        time              1   137  1.94  0.166    ""      0.014   
## 6 1        iu_group:time     1   137  0.816 0.368    ""      0.006
# results: 
# The effect of IU group in response to CS+ was significant [F(1,137) = 6.70, p = .011, pes = .047]
# the effect of IU group in response to CS- was also significant [F(1,137) = 14.43, p < .001, pes = .095]

ANOVA Extinction Saccade Amplitude

# transform wide format data into long format for mixed ANOVA 
df_long_ext_sacc_amplitude <- melt(df, id = c("id", "iu_group"), 
                                 measure.vars = c("e_ext_csp_sacc_amplitude", 
                                                  "e_ext_csm_sacc_amplitude",
                                                  "l_ext_csp_sacc_amplitude", 
                                                  "l_ext_csm_sacc_amplitude"))

# rename columns for easier interpretation
colnames(df_long_ext_sacc_amplitude) = c("id", "iu_group", "condition", "sacc_amplitude")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_ext_sacc_amplitude$stimulus <- 
  factor(ifelse(df_long_ext_sacc_amplitude$condition == "e_ext_csp_sacc_amplitude" | 
                  df_long_ext_sacc_amplitude$condition == "l_ext_csp_sacc_amplitude", 1, -1))

# create column to code extinction as early (1) and late (-1)
df_long_ext_sacc_amplitude$time <- 
  factor(ifelse(df_long_ext_sacc_amplitude$condition == "e_ext_csp_sacc_amplitude" | 
                  df_long_ext_sacc_amplitude$condition == "e_ext_csm_sacc_amplitude", 1, -1))

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) x 2 (Time: Early, Late) mixed ANOVA, 
# and obtain effect size (partial eta squared)
ext_sacc_amplitude_anova <- 
  anova_test(df_long_ext_sacc_amplitude, 
             sacc_amplitude ~ iu_group * stimulus * time + Error(id/(stimulus*time)),
             effect.size = "pes") 
## Warning: NA detected in rows: 116,181,301.
## Removing this rows before the analysis.
# obtain the mixed ANOVA results
get_anova_table(ext_sacc_amplitude_anova)
## ANOVA Table (type III tests)
## 
##                   Effect DFn DFd     F     p p<.05      pes
## 1               iu_group   1 134 3.170 0.077       0.023000
## 2               stimulus   1 134 0.740 0.391       0.005000
## 3                   time   1 134 0.275 0.601       0.002000
## 4      iu_group:stimulus   1 134 1.687 0.196       0.012000
## 5          iu_group:time   1 134 0.131 0.718       0.000977
## 6          stimulus:time   1 134 0.077 0.781       0.000577
## 7 iu_group:stimulus:time   1 134 0.609 0.437       0.005000
# results:
# IU: F(1,134) = 3.17, p = .077, eta2(partial) = .023
# Stimulus: F(1,134) = 0.74, p = .391, eta2(partial) = .005
# Time: F(1,134) = 0.28, p = .601, eta2(partial) = .002
# IU * Stimulus: F(1, 134) = 1.69, p = .196, eta2(partial) = .012
# IU * Time: F(1,134) = 0.13, p = .718, eta2(partial) < .001
# Stimulus * Time: F(1,134) = 0.08, p = .781, eta2(partial) < .001
# IU * Stimulus * Time: F(1,134) = .61, p = .437, eta2(partial) < .001

# therefore, there are no significant effects or interactions 
# on saccade amplitude throughout extinction 

# write to csv
write.csv((get_anova_table(ext_sacc_amplitude_anova)), 
          file = "tables/anovas/ext_sacc_amplitude_anova.csv")

Bar Graphs

Fixation Count

Acquisition

# obtain mean fix count for each group at each stimulus type and save as vector 
mean_acq_fix_count_high_iu_csp <- 
  mean(df$acq_csp_fix_count[df_long_acq_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS+
mean_acq_fix_count_low_iu_csp <- 
  mean(df$acq_csp_fix_count[df_long_acq_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS+
mean_acq_fix_count_high_iu_csm <- 
  mean(df$acq_csm_fix_count[df_long_acq_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS-
mean_acq_fix_count_low_iu_csm <- 
  mean(df$acq_csm_fix_count[df_long_acq_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS-

# combine into single variable called 
all_mean_acq_fix_count <- 
  c(mean_acq_fix_count_high_iu_csp, mean_acq_fix_count_low_iu_csp, 
    mean_acq_fix_count_high_iu_csm, mean_acq_fix_count_low_iu_csm)

# obtain SD fix count for each group at each stimulus type and save as vector 
sd_acq_fix_count_high_iu_csp <- 
  sd(df$acq_csp_fix_count[df_long_acq_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS+
sd_acq_fix_count_low_iu_csp <- 
  sd(df$acq_csp_fix_count[df_long_acq_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS+
sd_acq_fix_count_high_iu_csm <- 
  sd(df$acq_csm_fix_count[df_long_acq_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS-
sd_acq_fix_count_low_iu_csm <- 
  sd(df$acq_csm_fix_count[df_long_acq_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS-

# obtain SE:
se_acq_fix_count_high_iu_csp <- sd_acq_fix_count_high_iu_csp/sqrt(length(df$id))
se_acq_fix_count_low_iu_csp <- sd_acq_fix_count_low_iu_csp/sqrt(length(df$id))
se_acq_fix_count_high_iu_csm <- sd_acq_fix_count_high_iu_csm/sqrt(length(df$id))
se_acq_fix_count_low_iu_csm <- sd_acq_fix_count_low_iu_csm/sqrt(length(df$id))

# Combine all into single variable called all_se
all_se_acq_fix_count <- c(se_acq_fix_count_high_iu_csp, se_acq_fix_count_low_iu_csp,
                          se_acq_fix_count_high_iu_csm, se_acq_fix_count_low_iu_csm)

### Create new data frame for figures
# Which includes mean and SE for each condition
df_fig_acquisition_fix_count <- data.frame(all_mean_acq_fix_count, all_se_acq_fix_count)

### add labels 
# add two more variables to indicate IU group and stimulus type. 
# for IU group
df_fig_acquisition_fix_count$iu_group[1] <- "High IU"
df_fig_acquisition_fix_count$iu_group[2] <- "Low IU"
df_fig_acquisition_fix_count$iu_group[3] <- "High IU"
df_fig_acquisition_fix_count$iu_group[4] <- "Low IU"

# for stimulus
df_fig_acquisition_fix_count$stimulus[1] <- "CS+"
df_fig_acquisition_fix_count$stimulus[2] <- "CS+"
df_fig_acquisition_fix_count$stimulus[3] <- "CS-"
df_fig_acquisition_fix_count$stimulus[4] <- "CS-"

# and re-order levels of stimulus factor so that CS+ appears on left in the graph
df_fig_acquisition_fix_count$stimulus <- 
  factor(df_fig_acquisition_fix_count$stimulus,levels=c("CS+","CS-"))


### Create figure
 fig_acquisition_fix_count <- ggplot(df_fig_acquisition_fix_count, 
                                   aes(x = iu_group, y = all_mean_acq_fix_count, 
                                       fill = stimulus)) + 
   geom_bar(stat = "identity", position = position_dodge(.6), width = .5, alpha = .85) +
   scale_y_continuous(limits = c(0, 10), expand = c(0,0)) +
   theme_classic() + 
   theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
   theme(axis.text.y = element_text(size = 10),axis.ticks.y = element_line(size = 1), 
         axis.line.y = element_line(colour = "black")) +
   theme(axis.text.x = element_text(colour = "black"), axis.ticks.x = element_blank(),
         axis.line.x = element_line(colour = "black")) +
   theme(legend.position = "bottom", legend.title = element_blank()) +
   scale_fill_manual(values = c("#c45150", "#824372")) +
   ggtitle("Acquisition") +
   labs(y = "Mean Fixation Count", x = "Intolerance of Uncertainty") +
      geom_errorbar(aes(ymin = all_mean_acq_fix_count - all_se_acq_fix_count, 
                        ymax = all_mean_acq_fix_count + all_se_acq_fix_count), 
                 width = .15, position = position_dodge(.6), colour = "#090707", size = .3) 

print(fig_acquisition_fix_count) 

ggsave(filename = "graphs/bar_plots/acquisition_fix_count.png", 
       plot = fig_acquisition_fix_count, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Extinction

# obtain mean fix count for each group at each stimulus type and save as vector 
mean_e_ext_fix_count_high_iu_csp <- 
  mean(df$e_ext_csp_fix_count[df_long_ext_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS+ early
mean_e_ext_fix_count_low_iu_csp <- 
  mean(df$e_ext_csp_fix_count[df_long_ext_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS+ early
mean_e_ext_fix_count_high_iu_csm <- 
  mean(df$e_ext_csm_fix_count[df_long_ext_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS- early
mean_e_ext_fix_count_low_iu_csm <- 
  mean(df$e_ext_csm_fix_count[df_long_ext_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS- early
mean_l_ext_fix_count_high_iu_csp <- 
  mean(df$l_ext_csp_fix_count[df_long_ext_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS+ late
mean_l_ext_fix_count_low_iu_csp <- 
  mean(df$l_ext_csp_fix_count[df_long_ext_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS+ late
mean_l_ext_fix_count_high_iu_csm <- 
  mean(df$l_ext_csm_fix_count[df_long_ext_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS- late
mean_l_ext_fix_count_low_iu_csm <- 
  mean(df$l_ext_csm_fix_count[df_long_ext_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS- late

# combine into single variable called 
all_mean_ext_fix_count <- 
  c(mean_e_ext_fix_count_high_iu_csp, mean_e_ext_fix_count_low_iu_csp, 
    mean_e_ext_fix_count_high_iu_csm, mean_e_ext_fix_count_low_iu_csm,
    mean_l_ext_fix_count_high_iu_csp, mean_l_ext_fix_count_low_iu_csp, 
    mean_l_ext_fix_count_high_iu_csm, mean_l_ext_fix_count_low_iu_csm)

# obtain SD fix count for each group at each stimulus type and save as vector 
sd_e_ext_fix_count_high_iu_csp <- 
  sd(df$e_ext_csp_fix_count[df_long_ext_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS+ early
sd_e_ext_fix_count_low_iu_csp <- 
  sd(df$e_ext_csp_fix_count[df_long_ext_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS+ early
sd_e_ext_fix_count_high_iu_csm <- 
  sd(df$e_ext_csm_fix_count[df_long_ext_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS- early
sd_e_ext_fix_count_low_iu_csm <- 
  sd(df$e_ext_csm_fix_count[df_long_ext_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS- early
sd_l_ext_fix_count_high_iu_csp <- 
  sd(df$l_ext_csp_fix_count[df_long_ext_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS+ late
sd_l_ext_fix_count_low_iu_csp <- 
  sd(df$l_ext_csp_fix_count[df_long_ext_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS+ late
sd_l_ext_fix_count_high_iu_csm <- 
  sd(df$l_ext_csm_fix_count[df_long_ext_fix_count$iu_group == "1"], na.rm = TRUE) # high IU CS- late
sd_l_ext_fix_count_low_iu_csm <- 
  sd(df$l_ext_csm_fix_count[df_long_ext_fix_count$iu_group == "-1"], na.rm = TRUE) # low IU CS- late

# obtain SE:
se_e_ext_fix_count_high_iu_csp <- sd_e_ext_fix_count_high_iu_csp/sqrt(length(df$id))
se_e_ext_fix_count_low_iu_csp <- sd_e_ext_fix_count_low_iu_csp/sqrt(length(df$id))
se_e_ext_fix_count_high_iu_csm <- sd_e_ext_fix_count_high_iu_csm/sqrt(length(df$id))
se_e_ext_fix_count_low_iu_csm <- sd_e_ext_fix_count_low_iu_csm/sqrt(length(df$id))
se_l_ext_fix_count_high_iu_csp <- sd_l_ext_fix_count_high_iu_csp/sqrt(length(df$id))
se_l_ext_fix_count_low_iu_csp <- sd_l_ext_fix_count_low_iu_csp/sqrt(length(df$id))
se_l_ext_fix_count_high_iu_csm <- sd_l_ext_fix_count_high_iu_csm/sqrt(length(df$id))
se_l_ext_fix_count_low_iu_csm <- sd_l_ext_fix_count_low_iu_csm/sqrt(length(df$id))

# Combine all into single variable called all_se
all_se_ext_fix_count <- c(se_e_ext_fix_count_high_iu_csp, se_e_ext_fix_count_low_iu_csp,
                          se_e_ext_fix_count_high_iu_csm, se_e_ext_fix_count_low_iu_csm,
                          se_l_ext_fix_count_high_iu_csp, se_l_ext_fix_count_low_iu_csp, 
                          se_l_ext_fix_count_high_iu_csm, se_l_ext_fix_count_low_iu_csm)

### Create new data frame for figures
# Which includes mean and SE for each condition
df_fig_extinction_fix_count <- data.frame(all_mean_ext_fix_count, all_se_ext_fix_count)

### add labels 
# add two more variables to indicate IU group and stimulus type. 
# for IU group
df_fig_extinction_fix_count$iu_group[1] <- "High IU"
df_fig_extinction_fix_count$iu_group[2] <- "Low IU"
df_fig_extinction_fix_count$iu_group[3] <- "High IU"
df_fig_extinction_fix_count$iu_group[4] <- "Low IU"
df_fig_extinction_fix_count$iu_group[5] <- "High IU"
df_fig_extinction_fix_count$iu_group[6] <- "Low IU"
df_fig_extinction_fix_count$iu_group[7] <- "High IU"
df_fig_extinction_fix_count$iu_group[8] <- "Low IU"

# for stimulus
df_fig_extinction_fix_count$stimulus[1] <- "CS+"
df_fig_extinction_fix_count$stimulus[2] <- "CS+"
df_fig_extinction_fix_count$stimulus[3] <- "CS-"
df_fig_extinction_fix_count$stimulus[4] <- "CS-"
df_fig_extinction_fix_count$stimulus[5] <- "CS+"
df_fig_extinction_fix_count$stimulus[6] <- "CS+"
df_fig_extinction_fix_count$stimulus[7] <- "CS-"
df_fig_extinction_fix_count$stimulus[8] <- "CS-"

# and re-order levels of stimulus factor so that CS+ appears on left in the graph
df_fig_extinction_fix_count$stimulus <- 
  factor(df_fig_extinction_fix_count$stimulus,levels=c("CS+","CS-"))

# for early / late extinction
df_fig_extinction_fix_count$time[1] <- "Early"
df_fig_extinction_fix_count$time[2] <- "Early"
df_fig_extinction_fix_count$time[3] <- "Early"
df_fig_extinction_fix_count$time[4] <- "Early"
df_fig_extinction_fix_count$time[5] <- "Late"
df_fig_extinction_fix_count$time[6] <- "Late"
df_fig_extinction_fix_count$time[7] <- "Late"
df_fig_extinction_fix_count$time[8] <- "Late"

### create figure
 fig_extinction_fix_count <- ggplot(df_fig_extinction_fix_count, 
                                   aes(x = iu_group, y = all_mean_ext_fix_count, 
                                       fill = stimulus)) + 
   geom_bar(stat = "identity", position = position_dodge(.6), width = .5, alpha = .85) +
   scale_y_continuous(limits = c(0, 10), expand = c(0,0)) +
   facet_wrap(~ time) +
   theme_classic() + 
   theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
   theme(axis.text.y = element_text(size = 10),axis.ticks.y = element_line(size = 1), 
         axis.line.y = element_line(colour = "black")) +
   theme(axis.text.x = element_text(colour = "black"), axis.ticks.x = element_blank(),
         axis.line.x = element_line(colour = "black")) +
   theme(legend.position = "bottom", legend.title = element_blank()) +
   scale_fill_manual(values = c("#c45150", "#824372")) +
   ggtitle("Extinction") +
   labs(y = "Mean Fixation Count", x = "Intolerance of Uncertainty") +
      geom_errorbar(aes(ymin = all_mean_ext_fix_count - all_se_ext_fix_count, 
                        ymax = all_mean_ext_fix_count + all_se_ext_fix_count), 
                 width = .15, position = position_dodge(.6), colour = "#090707", size = .3) +
   theme(strip.background = element_blank()) +
  theme(strip.text = element_text(size = 12)) 

 # obtain and check figure 
print(fig_extinction_fix_count) 

# save figure to files
ggsave(filename = "graphs/bar_plots/extinction_fix_count.png", 
       plot = fig_extinction_fix_count, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Fixation Duration (Log Transformed)

Acquisition

# obtain mean fix duration for each group at each stimulus type and save as vector 
mean_acq_fix_duration_high_iu_csp_log <- 
  mean(df$acq_csp_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "1"], na.rm = TRUE) # high IU CS+
mean_acq_fix_duration_low_iu_csp_log <- 
  mean(df$acq_csp_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "-1"], na.rm = TRUE) # low IU CS+
mean_acq_fix_duration_high_iu_csm_log <- 
  mean(df$acq_csm_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "1"], na.rm = TRUE) # high IU CS-
mean_acq_fix_duration_low_iu_csm_log <- 
  mean(df$acq_csm_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "-1"], na.rm = TRUE) # low IU CS-

# combine into single variable called 
all_mean_acq_fix_duration_log <- 
  c(mean_acq_fix_duration_high_iu_csp_log, mean_acq_fix_duration_low_iu_csp_log, 
    mean_acq_fix_duration_high_iu_csm_log, mean_acq_fix_duration_low_iu_csm_log)

# obtain SD fix duration for each group at each stimulus type and save as vector 
sd_acq_fix_duration_high_iu_csp_log <- 
  sd(df$acq_csp_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "1"], na.rm = TRUE) # high IU CS+
sd_acq_fix_duration_low_iu_csp_log <- 
  sd(df$acq_csp_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "-1"], na.rm = TRUE) # low IU CS+
sd_acq_fix_duration_high_iu_csm_log <- 
  sd(df$acq_csm_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "1"], na.rm = TRUE) # high IU CS-
sd_acq_fix_duration_low_iu_csm_log <- 
  sd(df$acq_csm_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "-1"], na.rm = TRUE) # low IU CS-

# obtain SE:
se_acq_fix_duration_high_iu_csp_log <- sd_acq_fix_duration_high_iu_csp_log/sqrt(length(df$id))
se_acq_fix_duration_low_iu_csp_log <- sd_acq_fix_duration_low_iu_csp_log/sqrt(length(df$id))
se_acq_fix_duration_high_iu_csm_log <- sd_acq_fix_duration_high_iu_csm_log/sqrt(length(df$id))
se_acq_fix_duration_low_iu_csm_log <- sd_acq_fix_duration_low_iu_csm_log/sqrt(length(df$id))

# combine all into single variable called all_se
all_se_acq_fix_duration_log <- c(se_acq_fix_duration_high_iu_csp_log, se_acq_fix_duration_low_iu_csp_log,
                          se_acq_fix_duration_high_iu_csm_log, se_acq_fix_duration_low_iu_csm_log)

# create new data frame for figures, which includes mean and se 
# for each condition 
df_fig_acquisition_fix_duration_log <- data.frame(all_mean_acq_fix_duration_log, all_se_acq_fix_duration_log)

# add labels - add two more variables to indicate IU group and stimulus type. 
# for IU group
df_fig_acquisition_fix_duration_log$iu_group[1] <- "High IU"
df_fig_acquisition_fix_duration_log$iu_group[2] <- "Low IU"
df_fig_acquisition_fix_duration_log$iu_group[3] <- "High IU"
df_fig_acquisition_fix_duration_log$iu_group[4] <- "Low IU"

# for stimulus
df_fig_acquisition_fix_duration_log$stimulus[1] <- "CS+"
df_fig_acquisition_fix_duration_log$stimulus[2] <- "CS+"
df_fig_acquisition_fix_duration_log$stimulus[3] <- "CS-"
df_fig_acquisition_fix_duration_log$stimulus[4] <- "CS-"

# and re-order levels of stimulus factor so that CS+ appears on left in the graph
df_fig_acquisition_fix_duration_log$stimulus <- 
  factor(df_fig_acquisition_fix_duration_log$stimulus,levels = c("CS+","CS-"))

# create figure
 fig_acquisition_fix_duration_log <- ggplot(df_fig_acquisition_fix_duration_log, 
                                   aes(x = iu_group, y = all_mean_acq_fix_duration_log, 
                                       fill = stimulus)) + 
   geom_bar(stat = "identity", position = position_dodge(.6), width = .5, alpha = .85) +
   scale_y_continuous(limits = c(0, 7.2), expand = c(0,0)) +
   theme_classic() + 
   theme(text = element_text(family = "serif"), 
         plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
   theme(axis.text.y = element_text(size = 10),axis.ticks.y = element_line(size = 1), 
         axis.line.y = element_line(colour = "black")) +
   theme(axis.text.x = element_text(colour = "black"), axis.ticks.x = element_blank(),
         axis.line.x = element_line(colour = "black")) +
   theme(legend.position = "bottom", legend.title = element_blank()) +
   scale_fill_manual(values = c("#c45150", "#824372")) +
   ggtitle("Acquisition") +
   labs(y = "Mean Fixation Duration (ms) (Log Transformed)", x = "Intolerance of Uncertainty") +
      geom_errorbar(aes(ymin = all_mean_acq_fix_duration_log - all_se_acq_fix_duration_log, 
                        ymax = all_mean_acq_fix_duration_log + all_se_acq_fix_duration_log), 
                 width = .15, position = position_dodge(.6), colour = "#090707", size = .3) 

 # obtain and check figure 
print(fig_acquisition_fix_duration_log) 

# save figure to file 
ggsave(filename = "graphs/bar_plots/acquisition_fix_duration_log.png", 
       plot = fig_acquisition_fix_duration_log, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Extinction

# obtain mean fix duration for each group at each stimulus type and save as vector 
# high IU CS+ early
mean_e_ext_fix_duration_high_iu_csp_log <- 
  mean(df$e_ext_csp_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "1"], na.rm = TRUE) 

# low IU CS+ early
mean_e_ext_fix_duration_low_iu_csp_log <- 
  mean(df$e_ext_csp_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "-1"], na.rm = TRUE) 

# high IU CS- early
mean_e_ext_fix_duration_high_iu_csm_log <- 
  mean(df$e_ext_csm_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "1"], na.rm = TRUE) 

# low IU CS- early
mean_e_ext_fix_duration_low_iu_csm_log <- 
  mean(df$e_ext_csm_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "-1"], na.rm = TRUE) 

# high IU CS+ late
mean_l_ext_fix_duration_high_iu_csp_log <- 
  mean(df$l_ext_csp_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "1"], na.rm = TRUE) 

# low IU CS+ late
mean_l_ext_fix_duration_low_iu_csp_log <- 
  mean(df$l_ext_csp_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "-1"], na.rm = TRUE) 

# high IU CS- late
mean_l_ext_fix_duration_high_iu_csm_log <- 
  mean(df$l_ext_csm_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "1"], na.rm = TRUE) 

# low IU CS- late
mean_l_ext_fix_duration_low_iu_csm_log <- 
  mean(df$l_ext_csm_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "-1"], na.rm = TRUE) 

# combine into single variable called 
all_mean_ext_fix_duration_log <- 
  c(mean_e_ext_fix_duration_high_iu_csp_log, mean_e_ext_fix_duration_low_iu_csp_log, 
    mean_e_ext_fix_duration_high_iu_csm_log, mean_e_ext_fix_duration_low_iu_csm_log,
    mean_l_ext_fix_duration_high_iu_csp_log, mean_l_ext_fix_duration_low_iu_csp_log, 
    mean_l_ext_fix_duration_high_iu_csm_log, mean_l_ext_fix_duration_low_iu_csm_log)

# obtain SD fix duration for each group at each stimulus type and save as vector 
# high IU CS+ early
sd_e_ext_fix_duration_high_iu_csp_log <- 
  sd(df$e_ext_csp_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "1"], na.rm = TRUE) 

# low IU CS+ early
sd_e_ext_fix_duration_low_iu_csp_log <- 
  sd(df$e_ext_csp_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "-1"], na.rm = TRUE) 

# high IU CS- early
sd_e_ext_fix_duration_high_iu_csm_log <- 
  sd(df$e_ext_csm_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "1"], na.rm = TRUE) 

# low IU CS- early
sd_e_ext_fix_duration_low_iu_csm_log <- 
  sd(df$e_ext_csm_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "-1"], na.rm = TRUE) 

# high IU CS+ late
sd_l_ext_fix_duration_high_iu_csp_log <- 
  sd(df$l_ext_csp_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "1"], na.rm = TRUE) 

# low IU CS+ late
sd_l_ext_fix_duration_low_iu_csp_log <- 
  sd(df$l_ext_csp_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "-1"], na.rm = TRUE) 

# high IU CS- late
sd_l_ext_fix_duration_high_iu_csm_log <- 
  sd(df$l_ext_csm_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "1"], na.rm = TRUE) 

# low IU CS- late
sd_l_ext_fix_duration_low_iu_csm_log <- 
  sd(df$l_ext_csm_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "-1"], na.rm = TRUE) 

# obtain SE:
se_e_ext_fix_duration_high_iu_csp_log <- sd_e_ext_fix_duration_high_iu_csp_log/sqrt(length(df$id))
se_e_ext_fix_duration_low_iu_csp_log <- sd_e_ext_fix_duration_low_iu_csp_log/sqrt(length(df$id))
se_e_ext_fix_duration_high_iu_csm_log <- sd_e_ext_fix_duration_high_iu_csm_log/sqrt(length(df$id))
se_e_ext_fix_duration_low_iu_csm_log <- sd_e_ext_fix_duration_low_iu_csm_log/sqrt(length(df$id))
se_l_ext_fix_duration_high_iu_csp_log <- sd_l_ext_fix_duration_high_iu_csp_log/sqrt(length(df$id))
se_l_ext_fix_duration_low_iu_csp_log <- sd_l_ext_fix_duration_low_iu_csp_log/sqrt(length(df$id))
se_l_ext_fix_duration_high_iu_csm_log <- sd_l_ext_fix_duration_high_iu_csm_log/sqrt(length(df$id))
se_l_ext_fix_duration_low_iu_csm_log <- sd_l_ext_fix_duration_low_iu_csm_log/sqrt(length(df$id))

# combine all into single variable 
all_se_ext_fix_duration_log <- c(se_e_ext_fix_duration_high_iu_csp_log, se_e_ext_fix_duration_low_iu_csp_log,
                          se_e_ext_fix_duration_high_iu_csm_log, se_e_ext_fix_duration_low_iu_csm_log,
                          se_l_ext_fix_duration_high_iu_csp_log, se_l_ext_fix_duration_low_iu_csp_log, 
                          se_l_ext_fix_duration_high_iu_csm_log, se_l_ext_fix_duration_low_iu_csm_log)

# create new data frame for figures which includes mean and SE for each condition
df_fig_extinction_fix_duration_log <- data.frame(all_mean_ext_fix_duration_log, all_se_ext_fix_duration_log)

# add labels - add two more variables to indicate IU group, stimulus type and extinction time 
# for IU group
df_fig_extinction_fix_duration_log$iu_group[1] <- "High IU"
df_fig_extinction_fix_duration_log$iu_group[2] <- "Low IU"
df_fig_extinction_fix_duration_log$iu_group[3] <- "High IU"
df_fig_extinction_fix_duration_log$iu_group[4] <- "Low IU"
df_fig_extinction_fix_duration_log$iu_group[5] <- "High IU"
df_fig_extinction_fix_duration_log$iu_group[6] <- "Low IU"
df_fig_extinction_fix_duration_log$iu_group[7] <- "High IU"
df_fig_extinction_fix_duration_log$iu_group[8] <- "Low IU"

# for stimulus
df_fig_extinction_fix_duration_log$stimulus[1] <- "CS+"
df_fig_extinction_fix_duration_log$stimulus[2] <- "CS+"
df_fig_extinction_fix_duration_log$stimulus[3] <- "CS-"
df_fig_extinction_fix_duration_log$stimulus[4] <- "CS-"
df_fig_extinction_fix_duration_log$stimulus[5] <- "CS+"
df_fig_extinction_fix_duration_log$stimulus[6] <- "CS+"
df_fig_extinction_fix_duration_log$stimulus[7] <- "CS-"
df_fig_extinction_fix_duration_log$stimulus[8] <- "CS-"

# and re-order levels of stimulus factor so that CS+ appears on left in the graph
df_fig_extinction_fix_duration_log$stimulus <- 
  factor(df_fig_extinction_fix_duration_log$stimulus,levels=c("CS+","CS-"))

# for early / late extinction
df_fig_extinction_fix_duration_log$time[1] <- "Early"
df_fig_extinction_fix_duration_log$time[2] <- "Early"
df_fig_extinction_fix_duration_log$time[3] <- "Early"
df_fig_extinction_fix_duration_log$time[4] <- "Early"
df_fig_extinction_fix_duration_log$time[5] <- "Late"
df_fig_extinction_fix_duration_log$time[6] <- "Late"
df_fig_extinction_fix_duration_log$time[7] <- "Late"
df_fig_extinction_fix_duration_log$time[8] <- "Late"

# create figure
fig_extinction_fix_duration_log <- ggplot(df_fig_extinction_fix_duration_log, 
                                   aes(x = iu_group, y = all_mean_ext_fix_duration_log, 
                                       fill = stimulus)) + 
  geom_bar(stat = "identity", position = position_dodge(.6), width = .5, alpha = .85) +
  scale_y_continuous(limits = c(0, 7.2), expand = c(0,0)) +
  facet_wrap(~ time) +
  theme_classic() + 
  theme(text = element_text(family = "serif"), 
        plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  theme(axis.text.y = element_text(size = 10),axis.ticks.y = element_line(size = 1), 
        axis.line.y = element_line(colour = "black")) +
  theme(axis.text.x = element_text(colour = "black"), axis.ticks.x = element_blank(),
        axis.line.x = element_line(colour = "black")) +
  theme(legend.position = "bottom", legend.title = element_blank()) +
  scale_fill_manual(values = c("#c45150", "#824372")) +
  ggtitle("Extinction") +
  labs(y = "Mean Fixation Duration (ms) (Log Transformed)", x = "Intolerance of Uncertainty") +
  geom_errorbar(aes(ymin = all_mean_ext_fix_duration_log - all_se_ext_fix_duration_log, 
                    ymax = all_mean_ext_fix_duration_log + all_se_ext_fix_duration_log), 
                width = .15, position = position_dodge(.6), colour = "#090707", size = .3) +
  theme(strip.background = element_blank()) +
  theme(strip.text = element_text(size = 12))

# obtain and check figure 
print(fig_extinction_fix_duration_log) 

# save figure to files
ggsave(filename = "graphs/bar_plots/extinction_fix_duration_log.png", 
       plot = fig_extinction_fix_duration_log, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Saccade Amplitude

Acquisition

# obtain mean sacc amplitude for each group at each stimulus type and save as vector 
mean_acq_sacc_amplitude_high_iu_csp <- 
  mean(df$acq_csp_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "1"], na.rm = TRUE) # high IU CS+
mean_acq_sacc_amplitude_low_iu_csp <- 
  mean(df$acq_csp_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) # low IU CS+
mean_acq_sacc_amplitude_high_iu_csm <- 
  mean(df$acq_csm_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "1"], na.rm = TRUE) # high IU CS-
mean_acq_sacc_amplitude_low_iu_csm <- 
  mean(df$acq_csm_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) # low IU CS-

# combine into single variable called 
all_mean_acq_sacc_amplitude <- 
  c(mean_acq_sacc_amplitude_high_iu_csp, mean_acq_sacc_amplitude_low_iu_csp, 
    mean_acq_sacc_amplitude_high_iu_csm, mean_acq_sacc_amplitude_low_iu_csm)

# obtain SD sacc amplitude for each group at each stimulus type and save as vector 
sd_acq_sacc_amplitude_high_iu_csp <- 
  sd(df$acq_csp_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "1"], na.rm = TRUE) # high IU CS+
sd_acq_sacc_amplitude_low_iu_csp <- 
  sd(df$acq_csp_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) # low IU CS+
sd_acq_sacc_amplitude_high_iu_csm <- 
  sd(df$acq_csm_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "1"], na.rm = TRUE) # high IU CS-
sd_acq_sacc_amplitude_low_iu_csm <- 
  sd(df$acq_csm_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) # low IU CS-

# obtain SE:
se_acq_sacc_amplitude_high_iu_csp <- sd_acq_sacc_amplitude_high_iu_csp/sqrt(length(df$id))
se_acq_sacc_amplitude_low_iu_csp <- sd_acq_sacc_amplitude_low_iu_csp/sqrt(length(df$id))
se_acq_sacc_amplitude_high_iu_csm <- sd_acq_sacc_amplitude_high_iu_csm/sqrt(length(df$id))
se_acq_sacc_amplitude_low_iu_csm <- sd_acq_sacc_amplitude_low_iu_csm/sqrt(length(df$id))

# combine all into single variable called all_se
all_se_acq_sacc_amplitude <- c(se_acq_sacc_amplitude_high_iu_csp, se_acq_sacc_amplitude_low_iu_csp,
                             se_acq_sacc_amplitude_high_iu_csm, se_acq_sacc_amplitude_low_iu_csm)

# create new data frame for figures, which includes mean and se 
# for each condition 
df_fig_acquisition_sacc_amplitude <- data.frame(all_mean_acq_sacc_amplitude, all_se_acq_sacc_amplitude)

# add labels - add two more variables to indicate IU group and stimulus type. 
# for IU group
df_fig_acquisition_sacc_amplitude$iu_group[1] <- "High IU"
df_fig_acquisition_sacc_amplitude$iu_group[2] <- "Low IU"
df_fig_acquisition_sacc_amplitude$iu_group[3] <- "High IU"
df_fig_acquisition_sacc_amplitude$iu_group[4] <- "Low IU"

# for stimulus
df_fig_acquisition_sacc_amplitude$stimulus[1] <- "CS+"
df_fig_acquisition_sacc_amplitude$stimulus[2] <- "CS+"
df_fig_acquisition_sacc_amplitude$stimulus[3] <- "CS-"
df_fig_acquisition_sacc_amplitude$stimulus[4] <- "CS-"

# and re-order levels of stimulus factor so that CS+ appears on left in the graph
df_fig_acquisition_sacc_amplitude$stimulus <- 
  factor(df_fig_acquisition_sacc_amplitude$stimulus,levels=c("CS+","CS-"))

# create figure
fig_acquisition_sacc_amplitude <- ggplot(df_fig_acquisition_sacc_amplitude, 
                                       aes(x = iu_group, y = all_mean_acq_sacc_amplitude, 
                                           fill = stimulus)) + 
  geom_bar(stat = "identity", position = position_dodge(.6), width = .5, alpha = .85) +
  scale_y_continuous(limits = c(0, 4), expand = c(0,0)) +
  theme_classic() + 
  theme(text = element_text(family = "serif"), 
        plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  theme(axis.text.y = element_text(size = 10),axis.ticks.y = element_line(size = 1), 
        axis.line.y = element_line(colour = "black")) +
  theme(axis.text.x = element_text(colour = "black"), axis.ticks.x = element_blank(),
        axis.line.x = element_line(colour = "black")) +
  theme(legend.position = "bottom", legend.title = element_blank()) +
  scale_fill_manual(values = c("#c45150", "#824372")) +
  ggtitle("Acquisition") +
  labs(y = "Mean Saccade Amplitude (degrees/ms)", x = "Intolerance of Uncertainty") +
  geom_errorbar(aes(ymin = all_mean_acq_sacc_amplitude - all_se_acq_sacc_amplitude, 
                    ymax = all_mean_acq_sacc_amplitude + all_se_acq_sacc_amplitude), 
                width = .15, position = position_dodge(.6), colour = "#090707", size = .3) 

# obtain and check figure 
print(fig_acquisition_sacc_amplitude) 

# save figure to file 
ggsave(filename = "graphs/bar_plots/acquisition_sacc_amplitude.png", 
       plot = fig_acquisition_sacc_amplitude, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Extinction

# obtain mean sacc amplitude for each group at each stimulus type and save as vector 
# high IU CS+ early
mean_e_ext_sacc_amplitude_high_iu_csp <- 
  mean(df$e_ext_csp_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "1"], na.rm = TRUE) 

# low IU CS+ early
mean_e_ext_sacc_amplitude_low_iu_csp <- 
  mean(df$e_ext_csp_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) 

# high IU CS- early
mean_e_ext_sacc_amplitude_high_iu_csm <- 
  mean(df$e_ext_csm_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "1"], na.rm = TRUE) 

# low IU CS- early
mean_e_ext_sacc_amplitude_low_iu_csm <- 
  mean(df$e_ext_csm_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) 

# high IU CS+ late
mean_l_ext_sacc_amplitude_high_iu_csp <- 
  mean(df$l_ext_csp_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "1"], na.rm = TRUE) 

# low IU CS+ late
mean_l_ext_sacc_amplitude_low_iu_csp <- 
  mean(df$l_ext_csp_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) 

# high IU CS- late
mean_l_ext_sacc_amplitude_high_iu_csm <- 
  mean(df$l_ext_csm_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "1"], na.rm = TRUE) 

# low IU CS- late
mean_l_ext_sacc_amplitude_low_iu_csm <- 
  mean(df$l_ext_csm_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) 

# combine into single variable called 
all_mean_ext_sacc_amplitude <- 
  c(mean_e_ext_sacc_amplitude_high_iu_csp, mean_e_ext_sacc_amplitude_low_iu_csp, 
    mean_e_ext_sacc_amplitude_high_iu_csm, mean_e_ext_sacc_amplitude_low_iu_csm,
    mean_l_ext_sacc_amplitude_high_iu_csp, mean_l_ext_sacc_amplitude_low_iu_csp, 
    mean_l_ext_sacc_amplitude_high_iu_csm, mean_l_ext_sacc_amplitude_low_iu_csm)

# obtain SD sacc amplitude for each group at each stimulus type and save as vector 
# high IU CS+ early
sd_e_ext_sacc_amplitude_high_iu_csp <- 
  sd(df$e_ext_csp_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "1"], na.rm = TRUE) 

# low IU CS+ early
sd_e_ext_sacc_amplitude_low_iu_csp <- 
  sd(df$e_ext_csp_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) 

# high IU CS- early
sd_e_ext_sacc_amplitude_high_iu_csm <- 
  sd(df$e_ext_csm_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "1"], na.rm = TRUE) 

# low IU CS- early
sd_e_ext_sacc_amplitude_low_iu_csm <- 
  sd(df$e_ext_csm_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) 

# high IU CS+ late
sd_l_ext_sacc_amplitude_high_iu_csp <- 
  sd(df$l_ext_csp_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "1"], na.rm = TRUE) 

# low IU CS+ late
sd_l_ext_sacc_amplitude_low_iu_csp <- 
  sd(df$l_ext_csp_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) 

# high IU CS- late
sd_l_ext_sacc_amplitude_high_iu_csm <- 
  sd(df$l_ext_csm_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "1"], na.rm = TRUE) 

# low IU CS- late
sd_l_ext_sacc_amplitude_low_iu_csm <- 
  sd(df$l_ext_csm_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "-1"], na.rm = TRUE) 

# obtain SE:
se_e_ext_sacc_amplitude_high_iu_csp <- sd_e_ext_sacc_amplitude_high_iu_csp/sqrt(length(df$id))
se_e_ext_sacc_amplitude_low_iu_csp <- sd_e_ext_sacc_amplitude_low_iu_csp/sqrt(length(df$id))
se_e_ext_sacc_amplitude_high_iu_csm <- sd_e_ext_sacc_amplitude_high_iu_csm/sqrt(length(df$id))
se_e_ext_sacc_amplitude_low_iu_csm <- sd_e_ext_sacc_amplitude_low_iu_csm/sqrt(length(df$id))
se_l_ext_sacc_amplitude_high_iu_csp <- sd_l_ext_sacc_amplitude_high_iu_csp/sqrt(length(df$id))
se_l_ext_sacc_amplitude_low_iu_csp <- sd_l_ext_sacc_amplitude_low_iu_csp/sqrt(length(df$id))
se_l_ext_sacc_amplitude_high_iu_csm <- sd_l_ext_sacc_amplitude_high_iu_csm/sqrt(length(df$id))
se_l_ext_sacc_amplitude_low_iu_csm <- sd_l_ext_sacc_amplitude_low_iu_csm/sqrt(length(df$id))

# combine all into single variable 
all_se_ext_sacc_amplitude <- c(se_e_ext_sacc_amplitude_high_iu_csp, se_e_ext_sacc_amplitude_low_iu_csp,
                             se_e_ext_sacc_amplitude_high_iu_csm, se_e_ext_sacc_amplitude_low_iu_csm,
                             se_l_ext_sacc_amplitude_high_iu_csp, se_l_ext_sacc_amplitude_low_iu_csp, 
                             se_l_ext_sacc_amplitude_high_iu_csm, se_l_ext_sacc_amplitude_low_iu_csm)

# create new data frame for figures which includes mean and SE for each condition
df_fig_extinction_sacc_amplitude <- data.frame(all_mean_ext_sacc_amplitude, all_se_ext_sacc_amplitude)

# add labels - add two more variables to indicate IU group, stimulus type and extinction time 
# for IU group
df_fig_extinction_sacc_amplitude$iu_group[1] <- "High IU"
df_fig_extinction_sacc_amplitude$iu_group[2] <- "Low IU"
df_fig_extinction_sacc_amplitude$iu_group[3] <- "High IU"
df_fig_extinction_sacc_amplitude$iu_group[4] <- "Low IU"
df_fig_extinction_sacc_amplitude$iu_group[5] <- "High IU"
df_fig_extinction_sacc_amplitude$iu_group[6] <- "Low IU"
df_fig_extinction_sacc_amplitude$iu_group[7] <- "High IU"
df_fig_extinction_sacc_amplitude$iu_group[8] <- "Low IU"

# for stimulus
df_fig_extinction_sacc_amplitude$stimulus[1] <- "CS+"
df_fig_extinction_sacc_amplitude$stimulus[2] <- "CS+"
df_fig_extinction_sacc_amplitude$stimulus[3] <- "CS-"
df_fig_extinction_sacc_amplitude$stimulus[4] <- "CS-"
df_fig_extinction_sacc_amplitude$stimulus[5] <- "CS+"
df_fig_extinction_sacc_amplitude$stimulus[6] <- "CS+"
df_fig_extinction_sacc_amplitude$stimulus[7] <- "CS-"
df_fig_extinction_sacc_amplitude$stimulus[8] <- "CS-"

# and re-order levels of stimulus factor so that CS+ appears on left in the graph
df_fig_extinction_sacc_amplitude$stimulus <- 
  factor(df_fig_extinction_sacc_amplitude$stimulus,levels=c("CS+","CS-"))

# for early / late extinction
df_fig_extinction_sacc_amplitude$time[1] <- "Early"
df_fig_extinction_sacc_amplitude$time[2] <- "Early"
df_fig_extinction_sacc_amplitude$time[3] <- "Early"
df_fig_extinction_sacc_amplitude$time[4] <- "Early"
df_fig_extinction_sacc_amplitude$time[5] <- "Late"
df_fig_extinction_sacc_amplitude$time[6] <- "Late"
df_fig_extinction_sacc_amplitude$time[7] <- "Late"
df_fig_extinction_sacc_amplitude$time[8] <- "Late"

# create figure
fig_extinction_sacc_amplitude <- ggplot(df_fig_extinction_sacc_amplitude, 
                                      aes(x = iu_group, y = all_mean_ext_sacc_amplitude, 
                                          fill = stimulus)) + 
  geom_bar(stat = "identity", position = position_dodge(.6), width = .5, alpha = .85) +
  scale_y_continuous(limits = c(0, 4), expand = c(0,0)) +
  facet_wrap(~ time) +
  theme_classic() + 
  theme(text = element_text(family = "serif"), 
        plot.title = element_text(face = "bold", hjust = 0.5, size = 15)) +
  theme(axis.text.y = element_text(size = 10),axis.ticks.y = element_line(size = 1), 
        axis.line.y = element_line(colour = "black")) +
  theme(axis.text.x = element_text(colour = "black"), axis.ticks.x = element_blank(),
        axis.line.x = element_line(colour = "black")) +
  theme(legend.position = "bottom", legend.title = element_blank()) +
  scale_fill_manual(values = c("#c45150", "#824372")) +
  ggtitle("Extinction") +
  labs(y = "Mean Saccade Amplitude (degrees/ms)", x = "Intolerance of Uncertainty") +
  geom_errorbar(aes(ymin = all_mean_ext_sacc_amplitude - all_se_ext_sacc_amplitude, 
                    ymax = all_mean_ext_sacc_amplitude + all_se_ext_sacc_amplitude), 
                width = .15, position = position_dodge(.6), colour = "#090707", size = .3) +
  theme(strip.background = element_blank()) +
  theme(strip.text = element_text(size = 12))

# obtain and check figure 
print(fig_extinction_sacc_amplitude) 

# save figure to files
ggsave(filename = "graphs/bar_plots/extinction_sacc_amplitude.png", 
       plot = fig_extinction_sacc_amplitude, 
       width = 20, 
       height = 10, 
       dpi = 300, 
       units = "cm")

Combine Bar Plots

all_bar_plots <- grid.arrange(fig_acquisition_fix_count, fig_extinction_fix_count,
                              fig_acquisition_fix_duration_log, fig_extinction_fix_duration_log,
                              fig_acquisition_sacc_amplitude, fig_extinction_sacc_amplitude, 
                              ncol = 2)

# save figure to files
ggsave(filename = "graphs/bar_plots/all_bar_plots.png", 
       plot = all_bar_plots, 
       width = 40, 
       height = 30, 
       dpi = 300, 
       units = "cm")

ANCOVAs to test Specificity of IU over Trait Anxiety

ANCOVA Acquisition Fixation Count

# transform wide format data into long format for mixed ANCOVA 
df_long_acq_fix_count <- melt(df, id = c("id", "iu_group", "sticsa_total"), 
                                 measure.vars = c("acq_csp_fix_count", 
                                                  "acq_csm_fix_count"))

# rename columns for easier interpretation
colnames(df_long_acq_fix_count) = c("id", "iu_group", "sticsa_total", "condition", "fix_count")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_acq_fix_count$stimulus <- 
  factor(ifelse(df_long_acq_fix_count$condition == "acq_csp_fix_count", 1, -1))

# mean centre continuous covariate (STICSA)
# to apply mean centring, first obtain average sticsa scores for all participants, 
# and save as a variable
df_long_acq_fix_count$sticsa_total_avg <- mean(df_long_acq_fix_count$sticsa_total)

# next, subtract this average from all participants' sticsa scores,
# and save as a variable
df_long_acq_fix_count$sticsa_total_centred <- 
  df_long_acq_fix_count$sticsa_total - df_long_acq_fix_count$sticsa_total_avg 
# from this we have mean sticsa scores after centring 

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) mixed ANCOVA, 
# with mean-cenred STICSA as covariate
# and obtain effect size (partial eta squared)
acq_fix_count_ancova <- 
  anova_test(df_long_acq_fix_count, fix_count ~ iu_group * stimulus + Error(id/stimulus),
             covariate = sticsa_total_centred, effect.size = "pes") 

# obtain the mixed ANCOVA results
get_anova_table(acq_fix_count_ancova)
## ANOVA Table (type III tests)
## 
##                          Effect DFn DFd      F        p p<.05      pes
## 1          sticsa_total_centred   1 136  0.059 0.808000       0.000434
## 2                      iu_group   1 136  3.191 0.076000       0.023000
## 3                      stimulus   1 136 11.622 0.000858     * 0.079000
## 4 sticsa_total_centred:stimulus   1 136  1.845 0.177000       0.013000
## 5             iu_group:stimulus   1 136  1.230 0.269000       0.009000
# results:
# STICSA (centred): F(1,136) = 0.06, p = .808, eta2(partial) = < .001
# IU: F(1,136) = 3.19, p = .076, eta2(partial) = .023
# Stimulus: F(1,136) = 11.62, p < .001***, eta2(partial) = .079
# STICSA * Stimulus: F(1,136) = 1.85, p = .177, eta2(partial) = .013
# IU * Stimulus: F(1, 136) = 1.23, p = .269, eta2(partial) = .009

# therefore, after accounting for trait anxiety, IU no longer has a significant
# effect on fixation count in acquisition, but stimulus continues to have
# significant effect. IU*Stimulus interaction also remains non-significant, 
# even after controlling for trait anxiety.

# write to csv
write.csv((get_anova_table(acq_fix_count_ancova)), 
          file = "tables/ancovas/acq_fix_count_ancova.csv")

# as there was a significant main effect of stimulus, obtain estimated
# marginal means to be reported (temporarily from SPSS):

emmeans_acq_fix_count_ancova_csp <- 6.92
emmeans_acq_fix_count_ancova_csm <- 7.32

ANCOVA Acquisition Fixation Duration (Log Transformed)

# transform wide format data into long format for mixed ANCOVA 
df_long_acq_fix_duration_log <- melt(df, id = c("id", "iu_group", "sticsa_total"), 
                                 measure.vars = c("acq_csp_fix_duration_log", 
                                                  "acq_csm_fix_duration_log"))

# rename columns for easier interpretation
colnames(df_long_acq_fix_duration_log) = c("id", "iu_group", "sticsa_total", "condition", "fix_duration_log")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_acq_fix_duration_log$stimulus <- 
  factor(ifelse(df_long_acq_fix_duration_log$condition == "acq_csp_fix_duration_log", 1, -1))

# mean centre continuous covariate (STICSA)
# to apply mean centring, first obtain average sticsa scores for all participants, 
# and save as a variable
df_long_acq_fix_duration_log$sticsa_total_avg <- mean(df_long_acq_fix_duration_log$sticsa_total)

# next, subtract this average from all participants' sticsa scores,
# and save as a variable
df_long_acq_fix_duration_log$sticsa_total_centred <- 
  df_long_acq_fix_duration_log$sticsa_total - df_long_acq_fix_duration_log$sticsa_total_avg 
# from this we have mean sticsa scores after centring 

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) mixed ANCOVA, 
# with mean-cenred STICSA as covariate
# and obtain effect size (partial eta squared)
acq_fix_duration_ancova_log <- 
  anova_test(df_long_acq_fix_duration_log, fix_duration_log ~ iu_group * stimulus + Error(id/stimulus),
             covariate = sticsa_total_centred, effect.size = "pes") 

# obtain the mixed ANCOVA results
get_anova_table(acq_fix_duration_ancova_log)
## ANOVA Table (type III tests)
## 
##                          Effect DFn DFd     F     p p<.05   pes
## 1          sticsa_total_centred   1 136 0.268 0.606       0.002
## 2                      iu_group   1 136 3.890 0.051       0.028
## 3                      stimulus   1 136 2.935 0.089       0.021
## 4 sticsa_total_centred:stimulus   1 136 0.409 0.524       0.003
## 5             iu_group:stimulus   1 136 1.674 0.198       0.012
# results:
# STICSA (centred): F(1,136) = 0.27, p = .606, eta2(partial) = .002
# IU: F(1,136) = 3.89, p = .051, eta2(partial) = .028
# Stimulus: F(1,136) = 2.94, p = .089, eta2(partial) = .021
# STICSA * Stimulus: F(1,136) = 0.41, p = .524, eta2(partial) = .003
# IU * Stimulus: F(1, 136) = 1.67, p = .198, eta2(partial) = .012

# there are no significant effects or interactions on fixation duration in acquisition. 

# write to csv
write.csv((get_anova_table(acq_fix_duration_ancova_log)), 
          file = "tables/ancovas/acq_fix_duration_ancova_log.csv")

ANCOVA Acquisition Saccade Amplitude

# transform wide format data into long format for mixed ANCOVA 
df_long_acq_sacc_amplitude <- melt(df, id = c("id", "iu_group", "sticsa_total"), 
                                 measure.vars = c("acq_csp_sacc_amplitude", 
                                                  "acq_csm_sacc_amplitude"))

# rename columns for easier interpretation
colnames(df_long_acq_sacc_amplitude) = c("id", "iu_group", "sticsa_total", "condition", "sacc_amplitude")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_acq_sacc_amplitude$stimulus <- 
  factor(ifelse(df_long_acq_sacc_amplitude$condition == "acq_csp_sacc_amplitude", 1, -1))

# mean centre continuous covariate (STICSA)
# to apply mean centring, first obtain average sticsa scores for all participants, 
# and save as a variable
df_long_acq_sacc_amplitude$sticsa_total_avg <- mean(df_long_acq_sacc_amplitude$sticsa_total)

# next, subtract this average from all participants' sticsa scores,
# and save as a variable
df_long_acq_sacc_amplitude$sticsa_total_centred <- 
  df_long_acq_sacc_amplitude$sticsa_total - df_long_acq_sacc_amplitude$sticsa_total_avg 
# from this we have mean sticsa scores after centring 

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) mixed ANCOVA, 
# with mean-cenred STICSA as covariate
# and obtain effect size (partial eta squared)
acq_sacc_amplitude_ancova <- 
  anova_test(df_long_acq_sacc_amplitude, sacc_amplitude ~ iu_group * stimulus + Error(id/stimulus),
             covariate = sticsa_total_centred, effect.size = "pes") 
## Warning: NA detected in rows: 234,259.
## Removing this rows before the analysis.
# obtain the mixed ANCOVA results
get_anova_table(acq_sacc_amplitude_ancova)
## ANOVA Table (type III tests)
## 
##                          Effect DFn DFd     F     p p<.05       pes
## 1          sticsa_total_centred   1 134 0.007 0.935       0.0000503
## 2                      iu_group   1 134 2.128 0.147       0.0160000
## 3                      stimulus   1 134 0.943 0.333       0.0070000
## 4 sticsa_total_centred:stimulus   1 134 0.643 0.424       0.0050000
## 5             iu_group:stimulus   1 134 0.864 0.354       0.0060000
# results:
# STICSA (centred): F(1,134) = 0.01, p = .935, eta2(partial) < .001
# IU: F(1,134) = 2.13, p = .147, eta2(partial) = .016
# Stimulus: F(1,134) = 0.94, p = .333, eta2(partial) = .007
# STICSA * Stimulus: F(1,134) = 0.64, p = .424, eta2(partial) = .005
# IU * Stimulus: F(1, 134) = 0.86, p = .354, eta2(partial) = .006

# therefore, after accounting for trait anxiety, there continue not
# to be any significant effects of IU, stimulus, and interaction
# effects on saccade amplitude in acquisition. 

# write to csv
write.csv((get_anova_table(acq_sacc_amplitude_ancova)), 
          file = "tables/ancovas/acq_sacc_amplitude_ancova.csv")

ANCOVA Extinction Fixation Count

# transform wide format data into long format for mixed ANOVA 
df_long_ext_fix_count <- melt(df, id = c("id", "iu_group", "sticsa_total"), 
                                 measure.vars = c("e_ext_csp_fix_count", 
                                                  "e_ext_csm_fix_count",
                                                  "l_ext_csp_fix_count", 
                                                  "l_ext_csm_fix_count"))

# rename columns for easier interpretation
colnames(df_long_ext_fix_count) = c("id", "iu_group", "sticsa_total", "condition", "fix_count")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_ext_fix_count$stimulus <- 
  factor(ifelse(df_long_ext_fix_count$condition == "e_ext_csp_fix_count" | 
                  df_long_ext_fix_count$condition == "l_ext_csp_fix_count", 1, -1))

# create column to code extinction as early (1) and late (-1)
df_long_ext_fix_count$time <- 
  factor(ifelse(df_long_ext_fix_count$condition == "e_ext_csp_fix_count" | 
                  df_long_ext_fix_count$condition == "e_ext_csm_fix_count", 1, -1))

# mean centre continuous covariate (STICSA)
# to apply mean centring, first obtain average sticsa scores for all participants, 
# and save as a variable
df_long_ext_fix_count$sticsa_total_avg <- mean(df_long_ext_fix_count$sticsa_total)

# next, subtract this average from all participants' sticsa scores,
# and save as a variable
df_long_ext_fix_count$sticsa_total_centred <- 
  df_long_ext_fix_count$sticsa_total - df_long_ext_fix_count$sticsa_total_avg 
# from this we have mean sticsa scores after centring 

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) x 2 (Time: Early, Late) mixed ANOVA, 
# with mean-centred STICSA as covariate,
# and obtain effect size (partial eta squared)
ext_fix_count_ancova <- 
  anova_test(df_long_ext_fix_count, 
             fix_count ~ iu_group * stimulus * time + Error(id/(stimulus*time)),
             covariate = sticsa_total_centred, effect.size = "pes") 

# obtain the mixed ANCOVA results
get_anova_table(ext_fix_count_ancova)
## ANOVA Table (type III tests)
## 
##                                Effect DFn DFd        F     p p<.05        pes
## 1                sticsa_total_centred   1 136 0.433000 0.512       0.00300000
## 2                            iu_group   1 136 4.361000 0.039     * 0.03100000
## 3                            stimulus   1 136 4.209000 0.042     * 0.03000000
## 4                                time   1 136 5.692000 0.018     * 0.04000000
## 5       sticsa_total_centred:stimulus   1 136 1.098000 0.297       0.00800000
## 6                   iu_group:stimulus   1 136 4.560000 0.035     * 0.03200000
## 7           sticsa_total_centred:time   1 136 0.000429 0.984       0.00000316
## 8                       iu_group:time   1 136 3.489000 0.064       0.02500000
## 9                       stimulus:time   1 136 0.066000 0.797       0.00048800
## 10 sticsa_total_centred:stimulus:time   1 136 0.901000 0.344       0.00700000
## 11             iu_group:stimulus:time   1 136 0.044000 0.834       0.00032500
# results:
# STICSA (centred): F(1,136) = 0.43, p = .512, eta2(partial) = .003
# IU: F(1,136) = 4.36, p = .039*, eta2(partial) = .031
# Stimulus: F(1,136) = 4.21, p = .042*, eta2(partial) = .030
# Time: F(1,136) = 5.69, p = .018 *, eta2(partial) = .040
# STICSA * Stimulus: F(1,136) = 1.10, p = .297, eta2(partial) = .008
# IU * Stimulus: F(1, 136) = 4.56, p = .035*, eta2(partial) = .032
# STICSA* Time: F(1,136) = 0.00, p = .982, eta2(partial) < .001
# IU * Time: F(1,136) = 3.49, p = .064, eta2(partial) = . 025
# Stimulus * Time: F(1,136) = 0.07, p = .797, eta2(partial) < .001
# STICSA * Stimulus * Time: F(1,136) = 0.90, p = .344, eta2(partial) = .007
# IU * Stimulus * Time: F(1,136) = 0.04, p = .834, eta2(partial) < .001

# therefore, after accounting for trait anxiety, IU, Stimulus, and Time 
# continue to have a significant effect on fixation duration in acquisition. 
# there is no longer a significant interaction effect of IU*Time, 
# but there is now a significant interaction effect of IU*stimulus

# write to csv
write.csv((get_anova_table(ext_fix_count_ancova)), 
          file = "tables/ancovas/ext_fix_count_ancova.csv")

# as there was a significant IU*Stimulus interaction, conduct simple 
# main effects analysis:
## obtain effect of IU at each level of stimulus 
simple_effects_ext_fix_count_iu_ancova <- df_long_ext_fix_count %>%
  group_by(stimulus) %>%
  anova_test(dv = fix_count, wid = id, between = iu_group, within = time,
             covariate = sticsa_total_centred, effect.size = "pes") %>%
  get_anova_table() 

# get the output
simple_effects_ext_fix_count_iu_ancova
## # A tibble: 10 × 8
##    stimulus Effect                      DFn   DFd     F     p `p<.05`   pes
##  * <fct>    <chr>                     <dbl> <dbl> <dbl> <dbl> <chr>   <dbl>
##  1 -1       sticsa_total_centred          1   136 0.142 0.707 ""      0.001
##  2 -1       iu_group                      1   136 6.66  0.011 "*"     0.047
##  3 -1       time                          1   136 5.02  0.027 "*"     0.036
##  4 -1       sticsa_total_centred:time     1   136 0.369 0.545 ""      0.003
##  5 -1       iu_group:time                 1   136 2.25  0.136 ""      0.016
##  6 1        sticsa_total_centred          1   136 0.796 0.374 ""      0.006
##  7 1        iu_group                      1   136 2.16  0.143 ""      0.016
##  8 1        time                          1   136 2.86  0.093 ""      0.021
##  9 1        sticsa_total_centred:time     1   136 0.253 0.616 ""      0.002
## 10 1        iu_group:time                 1   136 2.40  0.124 ""      0.017
# results: 
# The effect of IU group on CS+ was not significant [F(1,136) = 2.17, p = .143, pes = .016]
# the effect of IU group on CS- was significant [F(1,136) = 6.66, p = .011, pes = .047]

# as there were a significant main effects of IU, stimulus and time, 
# as well as a sig IU-STimulus interaction, obtain estimated
# marginal means to be reported (temporarily from SPSS):

# IU
emmeans_ext_fix_count_ancova_high_iu <- 8.15
emmeans_ext_fix_count_ancova_low_iu <- 6.87

# stimulus
emmeans_ext_fix_count_ancova_csp <- 7.37
emmeans_ext_fix_count_ancova_csm <- 7.65

# time
emmeans_ext_fix_count_ancova_early <- 7.29
emmeans_ext_fix_count_ancova_late <- 7.72

# IU-Stimulus interaction
emmeans_ext_fix_count_ancova_high_iu_csp <- 7.83
emmeans_ext_fix_count_ancova_high_iu_csm <- 8.46
emmeans_ext_fix_count_ancova_low_iu_csp <- 6.90
emmeans_ext_fix_count_ancova_low_iu_csm <- 6.84

ANCOVA Extinction Fixation Duration (Log Transformed)

# transform wide format data into long format for mixed ANCOVA 
df_long_ext_fix_duration_log <- melt(df, id = c("id", "iu_group", "sticsa_total"), 
                                 measure.vars = c("e_ext_csp_fix_duration_log", 
                                                  "e_ext_csm_fix_duration_log",
                                                  "l_ext_csp_fix_duration_log", 
                                                  "l_ext_csm_fix_duration_log"))

# rename columns for easier interpretation
colnames(df_long_ext_fix_duration_log) = c("id", "iu_group", "sticsa_total", "condition", "fix_duration_log")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_ext_fix_duration_log$stimulus <- 
  factor(ifelse(df_long_ext_fix_duration_log$condition == "e_ext_csp_fix_duration_log" | 
                  df_long_ext_fix_duration_log$condition == "l_ext_csp_fix_duration_log", 1, -1))

# create column to code extinction as early (1) and late (-1)
df_long_ext_fix_duration_log$time <- 
  factor(ifelse(df_long_ext_fix_duration_log$condition == "e_ext_csp_fix_duration_log" | 
                  df_long_ext_fix_duration_log$condition == "e_ext_csm_fix_duration_log", 1, -1))

# mean centre continuous covariate (STICSA)
# to apply mean centring, first obtain average sticsa scores for all participants, 
# and save as a variable
df_long_ext_fix_duration_log$sticsa_total_avg <- mean(df_long_ext_fix_duration_log$sticsa_total)

# next, subtract this average from all participants' sticsa scores,
# and save as a variable
df_long_ext_fix_duration_log$sticsa_total_centred <- 
  df_long_ext_fix_duration_log$sticsa_total - df_long_ext_fix_duration_log$sticsa_total_avg 
# from this we have mean sticsa scores after centring 

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) x 2 (Time: Early, Late) mixed ANOVA, 
# with mean-centred STICSA as covariate,
# and obtain effect size (partial eta squared)
ext_fix_duration_ancova_log <- 
  anova_test(df_long_ext_fix_duration_log, 
             fix_duration_log ~ iu_group * stimulus * time + Error(id/(stimulus*time)),
             covariate = sticsa_total_centred, effect.size = "pes") 

# obtain the mixed ANCOVA results
get_anova_table(ext_fix_duration_ancova_log)
## ANOVA Table (type III tests)
## 
##                                Effect DFn DFd     F     p p<.05        pes
## 1                sticsa_total_centred   1 136 0.001 0.972       0.00000901
## 2                            iu_group   1 136 8.365 0.004     * 0.05800000
## 3                            stimulus   1 136 0.514 0.475       0.00400000
## 4                                time   1 136 4.358 0.039     * 0.03100000
## 5       sticsa_total_centred:stimulus   1 136 0.195 0.659       0.00100000
## 6                   iu_group:stimulus   1 136 5.357 0.022     * 0.03800000
## 7           sticsa_total_centred:time   1 136 0.329 0.567       0.00200000
## 8                       iu_group:time   1 136 0.501 0.480       0.00400000
## 9                       stimulus:time   1 136 0.174 0.677       0.00100000
## 10 sticsa_total_centred:stimulus:time   1 136 0.221 0.639       0.00200000
## 11             iu_group:stimulus:time   1 136 0.379 0.539       0.00300000
# results:
# STICSA (centred): F(1,136) = 0.01, p = .972, eta2(partial) < .001
# IU: F(1,136) = 8.37, p = .004**, eta2(partial) = .058
# Stimulus: F(1,136) = 0.51, p = .475, eta2(partial) = .004
# Time: F(1,136) = 4.36, p = .039*, eta2(partial) = .031
# STICSA * Stimulus: F(1,136) = 0.20, p = .659, eta2(partial) = .001
# IU * Stimulus: F(1, 136) = 5.36, p = .022*, eta2(partial) = .038
# STICSA* Time: F(1,136) = 0.33, p = .567, eta2(partial) = .002
# IU * Time: F(1,136) = 0.50, p = .480, eta2(partial) = . 004
# Stimulus * Time: F(1,136) = 0.17, p = .677, eta2(partial) = .001
# STICSA * Stimulus * Time: F(1,136) = 0.22, p = .639, eta2(partial) = .002
# IU * Stimulus * Time: F(1,136) = 0.34, p = .539, eta2(partial) = .003

# there were significant main effects of IU, time,
# and a significant IU-stimulus interaction on fixation duration in extinction,
# and no further main effects or interactions. 

# write to csv
write.csv((get_anova_table(ext_fix_duration_ancova_log)), 
          file = "tables/ancovas/ext_fix_duration_ancova_log.csv")

# as there was a significant IU*Stimulus interaction, conduct simple 
# main effects analysis:
## obtain effect of IU at each level of stimulus 
simple_effects_ext_fix_duration_iu_ancova <- df_long_ext_fix_duration_log %>%
  group_by(stimulus) %>%
  anova_test(dv = fix_duration_log, wid = id, between = iu_group, within = time,
             covariate = sticsa_total_centred, effect.size = "pes") %>%
  get_anova_table() 

# get the output
simple_effects_ext_fix_duration_iu_ancova
## # A tibble: 10 × 8
##    stimulus Effect                      DFn   DFd      F     p `p<.05`       pes
##  * <fct>    <chr>                     <dbl> <dbl>  <dbl> <dbl> <chr>       <dbl>
##  1 -1       sticsa_total_centred          1   136  0.008 0.928 ""      0.0000602
##  2 -1       iu_group                      1   136 11.2   0.001 "*"     0.076    
##  3 -1       time                          1   136  4.37  0.038 "*"     0.031    
##  4 -1       sticsa_total_centred:time     1   136  0.627 0.43  ""      0.005    
##  5 -1       iu_group:time                 1   136  0.061 0.806 ""      0.000445 
##  6 1        sticsa_total_centred          1   136  0.027 0.87  ""      0.000199 
##  7 1        iu_group                      1   136  4.70  0.032 "*"     0.033    
##  8 1        time                          1   136  1.93  0.167 ""      0.014    
##  9 1        sticsa_total_centred:time     1   136  0.036 0.849 ""      0.000266 
## 10 1        iu_group:time                 1   136  0.771 0.381 ""      0.006
# results: 
# The effect of IU group on CS+ was not significant [F(1,136) = 4.70, p = .032, pes = .033]
# the effect of IU group on CS- was significant [F(1,136) = 11,19, p = .001, pes = .076]

# as there were a significant main effects of IU, time and iu-stimulus interaction, 
# obtain estimated marginal means to be reported (temporarily from SPSS):
##### NOTE: have not done log-transformed fixation duration on SPSS yet, and
# so do not have estimated marginal means for fixation duration (log-transformed)

# IU
emmeans_ext_fix_duration_ancova_high_iu <- 0
emmeans_ext_fix_duration_ancova_low_iu <- 0

# time
emmeans_ext_fix_duration_ancova_early <- 0
emmeans_ext_fix_duration_ancova_late <- 0

# IU-Stimulus interaction
emmeans_ext_fix_duration_ancova_high_iu_csp <- 0
emmeans_ext_fix_duration_ancova_high_iu_csm <- 0
emmeans_ext_fix_duration_ancova_low_iu_csp <- 0
emmeans_ext_fix_duration_ancova_low_iu_csm <- 0

ANCOVA Extinction Saccade Amplitude

# transform wide format data into long format for mixed ANCOVA 
df_long_ext_sacc_amplitude <- melt(df, id = c("id", "iu_group", "sticsa_total"), 
                                 measure.vars = c("e_ext_csp_sacc_amplitude", 
                                                  "e_ext_csm_sacc_amplitude",
                                                  "l_ext_csp_sacc_amplitude", 
                                                  "l_ext_csm_sacc_amplitude"))

# rename columns for easier interpretation
colnames(df_long_ext_sacc_amplitude) = c("id", "iu_group", "sticsa_total", "condition", "sacc_amplitude")

# create column to code stimulus as CS+ (1) and CS- (-1)
df_long_ext_sacc_amplitude$stimulus <- 
  factor(ifelse(df_long_ext_sacc_amplitude$condition == "e_ext_csp_sacc_amplitude" | 
                  df_long_ext_sacc_amplitude$condition == "l_ext_csp_sacc_amplitude", 1, -1))

# create column to code extinction as early (1) and late (-1)
df_long_ext_sacc_amplitude$time <- 
  factor(ifelse(df_long_ext_sacc_amplitude$condition == "e_ext_csp_sacc_amplitude" | 
                  df_long_ext_sacc_amplitude$condition == "e_ext_csm_sacc_amplitude", 1, -1))

# mean centre continuous covariate (STICSA)
# to apply mean centring, first obtain average sticsa scores for all participants, 
# and save as a variable
df_long_ext_sacc_amplitude$sticsa_total_avg <- mean(df_long_ext_sacc_amplitude$sticsa_total)

# next, subtract this average from all participants' sticsa scores,
# and save as a variable
df_long_ext_sacc_amplitude$sticsa_total_centred <- 
  df_long_ext_sacc_amplitude$sticsa_total - df_long_ext_sacc_amplitude$sticsa_total_avg 
# from this we have mean sticsa scores after centring 

# compute 2(IU: High & Low) x 2 (Stimulus: CS+, CS-) x 2 (Time: Early, Late) mixed ANOVA, 
# with mean-centred STICSA as covariate,
# and obtain effect size (partial eta squared)
ext_sacc_amplitude_ancova <- 
  anova_test(df_long_ext_sacc_amplitude, 
             sacc_amplitude ~ iu_group * stimulus * time + Error(id/(stimulus*time)),
             covariate = sticsa_total_centred, effect.size = "pes") 
## Warning: NA detected in rows: 116,181,301.
## Removing this rows before the analysis.
# obtain the mixed ANCOVA results
get_anova_table(ext_sacc_amplitude_ancova)
## ANOVA Table (type III tests)
## 
##                                Effect DFn DFd     F     p p<.05      pes
## 1                sticsa_total_centred   1 133 1.134 0.289       0.008000
## 2                            iu_group   1 133 1.025 0.313       0.008000
## 3                            stimulus   1 133 0.754 0.387       0.006000
## 4                                time   1 133 0.255 0.615       0.002000
## 5       sticsa_total_centred:stimulus   1 133 0.370 0.544       0.003000
## 6                   iu_group:stimulus   1 133 2.035 0.156       0.015000
## 7           sticsa_total_centred:time   1 133 1.359 0.246       0.010000
## 8                       iu_group:time   1 133 0.803 0.372       0.006000
## 9                       stimulus:time   1 133 0.071 0.790       0.000533
## 10 sticsa_total_centred:stimulus:time   1 133 0.421 0.517       0.003000
## 11             iu_group:stimulus:time   1 133 0.997 0.320       0.007000
# results:
# STICSA (centred): F(1,133) = 1.13, p = .289, eta2(partial) = .008
# IU: F(1,133) = 1.03, p = .313, eta2(partial) = .008
# Stimulus: F(1,133) = 0.75, p = .387, eta2(partial) = .006
# Time: F(1,133) = 0.26, p = .615, eta2(partial) = .002
# STICSA * Stimulus: F(1,133) = 0.37, p = .544, eta2(partial) = .003
# IU * Stimulus: F(1, 133) = 2.04, p = .156, eta2(partial) = .015
# STICSA* Time: F(1,133) = 1.36, p = .246, eta2(partial) = .010
# IU * Time: F(1,133) = 0.80, p = .372, eta2(partial) = . 006
# Stimulus * Time: F(1,133) = 0.07, p = .790, eta2(partial) = .001
# STICSA * Stimulus * Time: F(1,133) = 0.42, p = .517, eta2(partial) = .003
# IU * Stimulus * Time: F(1,133) = 0.10, p = .320, eta2(partial) = .007

# therefore, even after accounting for trait anxiety, there continue
# to be no significant effects or interactions on saccade amplitude
# in extinction

# write to csv
write.csv((get_anova_table(ext_sacc_amplitude_ancova)), 
          file = "tables/ancovas/ext_sacc_amplitude_ancova.csv")

Assumption Checks

# assumptions of mixed ANOVA:
## categorical IVs, interval/ratio DVs
## residuals approximate normal distribution (for each level of each IV)
## homogeneity of variances
## homoscedasticity (plot standardised residuals against predicted values)
## sphericity (not applicable in this case, as no within-subjects factors with > 3 levels)
## homogeneity of variance-covariance matrices (this is assumed as sample size for
##### each group and testing session is roughly equal)

# additional assumptions of ANCOVA:
## independence of covariate and IVs
## homogeneity of regression slopes
## linearity between covariate and outcome variable(s) at each level of grouping
#### variable(s). check by creating grouped scatterplot
## outcome variable(s) should be approximately normally distributed
## no sigifnicant outliers in the groups

Linearity Between Covariate and Outcome Variables

## this is at each level of grouping variable. 
# check by computing grouped scatterplot of covariate and outcome variable 

################# linearity 
# linearity between covariate and outcome variable at each level of grouping variable
# (check using grouped scatterplot of covariate and outcome variable)

### iu group
scatterplot_acq_fix_count_sticsa_centred_by_iu <- ggplot(df_long_acq_fix_count, 
                                       aes(x = sticsa_total_centred, y = fix_count, 
                                           colour = iu_group)) +
  geom_point() +
  geom_jitter(width = .5, alpha = .30, size = 2.5) +
  geom_smooth(method = lm, se = FALSE) +
  labs(title = "Plot of the Relationship Between Trait Anxiety (Covariate) and 
  Fixation Count in Acquisition by IU Group",
       x = "STICSA Total (Centred)", 
       y = "Acquisition Fixation Count") +
  theme_classic() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  theme(text = element_text(family = "serif")) +
  guides(colour = guide_legend(reverse = TRUE)) +
  scale_colour_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
    labs(colour = "IU Group") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) 

print(scatterplot_acq_fix_count_sticsa_centred_by_iu)
## `geom_smooth()` using formula 'y ~ x'

# there appears to be a linear relationship between the covariate (sticsa total centred)
# and outcome variable (fixation count acquisition) at both levels of the grouping variable (iu group)
# however, this plot also demonstrates likely heterogeneity of regression slopes, as there is possibly an
# interaction between the outcome and covariate, and the regression lines do not appear
# parallel

####### extinctino iu
### iu group
scatterplot_ext_fix_count_sticsa_centred_by_iu <- ggplot(df_long_ext_fix_count, 
                                       aes(x = sticsa_total_centred, y = fix_count, 
                                           colour = iu_group)) +
  geom_point() +
  geom_jitter(width = .5, alpha = .30, size = 2.5) +
  geom_smooth(method = lm, se = FALSE) +
  labs(title = "Plot of the Relationship Between Trait Anxiety (Covariate) and 
  Fixation Count in Extinction by IU Group and Time",
       x = "STICSA Total (Centred)", 
       y = "Extinction Fixation Count") +
  theme_classic() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  theme(text = element_text(family = "serif")) +
  guides(colour = guide_legend(reverse = TRUE)) +
  scale_colour_manual(values = c("#c45150", "#824372"), labels = c("Low IU", "High IU")) +
    labs(colour = "IU Group") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold")) +
  facet_wrap(~ time, labeller = "label_both") 

print(scatterplot_ext_fix_count_sticsa_centred_by_iu)
## `geom_smooth()` using formula 'y ~ x'

### stimulus
scatterplot_acq_fix_count_sticsa_centred_by_stimulus <- ggplot(df_long_acq_fix_count, 
                                       aes(x = sticsa_total_centred, y = fix_count, 
                                           colour = stimulus)) +
  geom_point() +
  geom_jitter(width = .5, alpha = .30, size = 2.5) +
  geom_smooth(method = lm, se = FALSE) +
  labs(title = "Plot of the Relationship Between Trait Anxiety (Covariate) and 
  Fixation Count in Extinction by Stimulus",
       x = "STICSA Total (Centred)", 
       y = "Acquisition Fixation Count") +
  theme_classic() +
  theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12)) +
  theme(text = element_text(family = "serif")) +
  guides(colour = guide_legend(reverse = TRUE)) +
  scale_colour_manual(values = c("#DCECC9", "#076B68"), labels = c("CS-", "CS+")) +
    labs(colour = "Stimulus") +
  theme(legend.position = "bottom", legend.title = element_text(face = "bold"))

print(scatterplot_acq_fix_count_sticsa_centred_by_stimulus)
## `geom_smooth()` using formula 'y ~ x'

# there appears to be a linear relationship between the covariate (sticsa total centred)
# and outcome variable (fixation count acquisition) at both levels of the grouping variable (stimulus)
# this plot also demonstrates homogeneity of regression slopes, as there is no 
# interaction between the outcome and covariate, and the regression lines appear roughly
# parallel

# this also checks homogeneity of regression slopes:
# slopes of regression lines should be same for each group (no interaction between
# covariate and outcome) i.e. the lines should be parallel

## reporting example: There was a linear relationship between pre-test and post-test anxiety score for each training group, as assessed by visual inspection of a scatter plot.

Normality

## acquisition fix count
df_long_acq_fix_count %>%
  group_by(iu_group, stimulus) %>%
  shapiro_test(fix_count)
## # A tibble: 4 × 5
##   iu_group stimulus variable  statistic         p
##   <fct>    <fct>    <chr>         <dbl>     <dbl>
## 1 -1       -1       fix_count     0.929 0.000611 
## 2 -1       1        fix_count     0.934 0.00108  
## 3 1        -1       fix_count     0.962 0.0364   
## 4 1        1        fix_count     0.895 0.0000312
## extinction fix count
df_long_ext_fix_count %>%
  group_by(iu_group, stimulus, time) %>%
  shapiro_test(fix_count)
## # A tibble: 8 × 6
##   iu_group stimulus time  variable  statistic          p
##   <fct>    <fct>    <fct> <chr>         <dbl>      <dbl>
## 1 -1       -1       -1    fix_count     0.961 0.0263    
## 2 -1       -1       1     fix_count     0.904 0.0000488 
## 3 -1       1        -1    fix_count     0.977 0.228     
## 4 -1       1        1     fix_count     0.881 0.00000681
## 5 1        -1       -1    fix_count     0.929 0.000810  
## 6 1        -1       1     fix_count     0.981 0.391     
## 7 1        1        -1    fix_count     0.945 0.00457   
## 8 1        1        1     fix_count     0.931 0.000995
# acquisition fix duration log
df_long_acq_fix_duration_log %>%
  group_by(iu_group, stimulus) %>%
  shapiro_test(fix_duration_log)
## # A tibble: 4 × 5
##   iu_group stimulus variable         statistic      p
##   <fct>    <fct>    <chr>                <dbl>  <dbl>
## 1 -1       -1       fix_duration_log     0.970 0.0814
## 2 -1       1        fix_duration_log     0.964 0.0398
## 3 1        -1       fix_duration_log     0.981 0.385 
## 4 1        1        fix_duration_log     0.981 0.408
## extinction fix count
df_long_ext_fix_duration_log %>%
  group_by(iu_group, stimulus, time) %>%
  shapiro_test(fix_duration_log)
## # A tibble: 8 × 6
##   iu_group stimulus time  variable         statistic       p
##   <fct>    <fct>    <fct> <chr>                <dbl>   <dbl>
## 1 -1       -1       -1    fix_duration_log     0.945 0.00364
## 2 -1       -1       1     fix_duration_log     0.974 0.143  
## 3 -1       1        -1    fix_duration_log     0.972 0.112  
## 4 -1       1        1     fix_duration_log     0.970 0.0913 
## 5 1        -1       -1    fix_duration_log     0.973 0.146  
## 6 1        -1       1     fix_duration_log     0.959 0.0242 
## 7 1        1        -1    fix_duration_log     0.984 0.523  
## 8 1        1        1     fix_duration_log     0.983 0.460
## acquisition sacc amplitude
df_long_acq_sacc_amplitude %>%
  group_by(iu_group, stimulus) %>%
  shapiro_test(sacc_amplitude)
## # A tibble: 4 × 5
##   iu_group stimulus variable       statistic        p
##   <fct>    <fct>    <chr>              <dbl>    <dbl>
## 1 -1       -1       sacc_amplitude     0.940 0.00227 
## 2 -1       1        sacc_amplitude     0.954 0.0111  
## 3 1        -1       sacc_amplitude     0.913 0.000176
## 4 1        1        sacc_amplitude     0.918 0.000275
## extinction sacc amplitude
df_long_ext_sacc_amplitude %>%
  group_by(iu_group, stimulus, time) %>%
  shapiro_test(sacc_amplitude)
## # A tibble: 8 × 6
##   iu_group stimulus time  variable       statistic            p
##   <fct>    <fct>    <fct> <chr>              <dbl>        <dbl>
## 1 -1       -1       -1    sacc_amplitude     0.849 0.000000535 
## 2 -1       -1       1     sacc_amplitude     0.889 0.0000125   
## 3 -1       1        -1    sacc_amplitude     0.821 0.0000000925
## 4 -1       1        1     sacc_amplitude     0.880 0.00000688  
## 5 1        -1       -1    sacc_amplitude     0.930 0.000946    
## 6 1        -1       1     sacc_amplitude     0.902 0.0000659   
## 7 1        1        -1    sacc_amplitude     0.925 0.000514    
## 8 1        1        1     sacc_amplitude     0.926 0.000578
## however, as sample size is greater than 50, normal QQ plot is preferred because
# at larger sample sizes, Shapiro-Wilk becomes very sensitive to even a minor 
# deviation from normality. 
# (q-q plot draw correlation btw data and the normal distribution)

ggqqplot(df_long_acq_fix_count, "fix_count", ggtheme = theme_classic()) +
           facet_grid(stimulus ~ iu_group, labeller = "label_both")

ggqqplot(df_long_ext_fix_count, "fix_count", ggtheme = theme_classic()) +
           facet_grid(stimulus + time ~ iu_group, labeller = "label_both")

### reporting example:
# The score were normally distributed (p > 0.05) for each cell, as assessed by Shapiro-Wilk’s test of normality.
# All the points fall approximately along the reference line, for each cell. So we can assume normality of the data.

Normality of Residuals

# first compute model using lm(). covariate goes first.f
ancova_lm_model_acq_fix_count <- 
  lm(fix_count ~ sticsa_total_centred + iu_group * stimulus, df_long_acq_fix_count)

# augment data by adding fitted values and residuals
model_metrics <- augment(ancova_lm_model_acq_fix_count) 

# assess normality of residuals using shapiro wilk
shapiro_test(model_metrics$.resid) # not normally distributed
## # A tibble: 1 × 3
##   variable             statistic       p.value
##   <chr>                    <dbl>         <dbl>
## 1 model_metrics$.resid     0.940 0.00000000340
# reporting:
# The Shapiro Wilk test was not significant (p > 0.05), so we can assume normality of residuals

###### also outliers:
# examine standardised/studentised residual. observations whose standardised residuals
# are greater than 3 in absolute value are possible outliers. 
model_metrics %>%
  filter(abs(.std.resid) > 3) %>%
  as.data.frame()
##   fix_count sticsa_total_centred iu_group stimulus  .fitted   .resid       .hat
## 1  23.16667             27.46043        1        1 7.697435 15.46923 0.04169528
## 2  20.16667             11.46043        1        1 7.564182 12.60248 0.01706000
## 3  18.33333             27.46043        1       -1 8.159690 10.17364 0.04169528
##     .sigma    .cooksd .std.resid
## 1 3.286205 0.18613043   4.624897
## 2 3.335139 0.04804383   3.720299
## 3 3.364537 0.08050687   3.041654
# reporting example: There were no outliers in the data, as assessed by no cases with standardized residuals greater than 3 in absolute value.

Homogeneity of Variance

# this will be done using levene's test
# however, in large samples, levene's test can be sig even when group variances
# are not very different.

## acquisition fix count
df_long_acq_fix_count %>%
  group_by(stimulus) %>%
  levene_test(fix_count ~ iu_group)
## # A tibble: 2 × 5
##   stimulus   df1   df2 statistic     p
##   <fct>    <int> <int>     <dbl> <dbl>
## 1 -1           1   137    0.477  0.491
## 2 1            1   137    0.0415 0.839
## extinction fix count
df_long_ext_fix_count %>%
  group_by(stimulus, time) %>%
  levene_test(fix_count ~ iu_group)
## # A tibble: 4 × 6
##   stimulus time    df1   df2 statistic     p
##   <fct>    <fct> <int> <int>     <dbl> <dbl>
## 1 -1       -1        1   137     1.45  0.231
## 2 -1       1         1   137     0.181 0.671
## 3 1        -1        1   137     0.264 0.608
## 4 1        1         1   137     1.86  0.174
# acquisition fix duration log
df_long_acq_fix_duration_log %>%
  group_by(stimulus) %>%
  levene_test(fix_duration_log ~ iu_group)
## # A tibble: 2 × 5
##   stimulus   df1   df2 statistic     p
##   <fct>    <int> <int>     <dbl> <dbl>
## 1 -1           1   137     2.04  0.155
## 2 1            1   137     0.753 0.387
## extinction fix count
df_long_ext_fix_duration_log %>%
  group_by(stimulus, time) %>%
  levene_test(fix_duration_log ~ iu_group)
## # A tibble: 4 × 6
##   stimulus time    df1   df2 statistic       p
##   <fct>    <fct> <int> <int>     <dbl>   <dbl>
## 1 -1       -1        1   137      8.18 0.00490
## 2 -1       1         1   137      7.74 0.00616
## 3 1        -1        1   137      2.78 0.0977 
## 4 1        1         1   137      7.14 0.00843
## acquisition sacc amplitude
df_long_acq_sacc_amplitude %>%
  group_by(stimulus) %>%
  levene_test(sacc_amplitude ~ iu_group)
## # A tibble: 2 × 5
##   stimulus   df1   df2 statistic      p
##   <fct>    <int> <int>     <dbl>  <dbl>
## 1 -1           1   135      1.03 0.311 
## 2 1            1   137      3.42 0.0665
## extinction sacc amplitude
df_long_ext_sacc_amplitude %>%
  group_by(stimulus, time) %>%
  levene_test(sacc_amplitude ~ iu_group)
## # A tibble: 4 × 6
##   stimulus time    df1   df2 statistic     p
##   <fct>    <fct> <int> <int>     <dbl> <dbl>
## 1 -1       -1        1   137    0.364  0.547
## 2 -1       1         1   136    1.72   0.191
## 3 1        -1        1   136    0.0230 0.880
## 4 1        1         1   136    0.0324 0.857
### reproting: For IUS total, the variances were not similar for the 
# High and Low IU groups, F(1,137) = 6.75, p = .010
# (i.e. spread of IU scores is different between groups)

# There was homogeneity of variances, as assessed by Levene’s test (p > 0.05).

Homoscedasticity

# aka homogeneity of residuals variance for all groups. the residuals are 
# assumed to have a constant variance (homoscedasticity)

Homogeneity of Variance-Covariance Matrices

# this tests whether covariance matrices are equal across cells formed by
# between-subjects factor (IU)

# use Box's M (however, this is highly sensitive, so unless p < .001 and sample
# sizes are unequal, can ignore it)

box_m(df_long_acq_fix_count[, "fix_count", drop = FALSE], df_long_acq_fix_count$iu_group)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1     0.224   0.636         1 Box's M-test for Homogeneity of Covariance Matric…
box_m(df_long_ext_fix_count[, "fix_count", drop = FALSE], df_long_ext_fix_count$iu_group)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1     0.753   0.385         1 Box's M-test for Homogeneity of Covariance Matric…
box_m(df_long_acq_fix_duration_log[, "fix_duration_log", drop = FALSE], df_long_acq_fix_duration_log$iu_group)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1     0.358   0.550         1 Box's M-test for Homogeneity of Covariance Matric…
box_m(df_long_ext_fix_duration_log[, "fix_duration_log", drop = FALSE], df_long_ext_fix_duration_log$iu_group)
## # A tibble: 1 × 4
##   statistic   p.value parameter method                                          
##       <dbl>     <dbl>     <dbl> <chr>                                           
## 1      16.7 0.0000435         1 Box's M-test for Homogeneity of Covariance Matr…
box_m(df_long_acq_sacc_amplitude[, "sacc_amplitude", drop = FALSE], df_long_acq_sacc_amplitude$iu_group)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1        NA      NA         1 Box's M-test for Homogeneity of Covariance Matric…
box_m(df_long_ext_sacc_amplitude[, "sacc_amplitude", drop = FALSE], df_long_ext_sacc_amplitude$iu_group)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1        NA      NA         1 Box's M-test for Homogeneity of Covariance Matric…
## reporting example: There was homogeneity of covariances, as assessed by Box’s test of equality of covariance matrices (p > 0.001).

Independence of Covariate and IVs

Fixation Count

Acquisition

# sticsa and iu group
t_test_independence_sticsa_iu_group_acq_fix_count <- 
  t.test(
    df_long_acq_fix_count[df_long_acq_fix_count$iu_group == "1", "sticsa_total_centred"],
    df_long_acq_fix_count[df_long_acq_fix_count$iu_group == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_iu_group_acq_fix_count
## 
##  Two Sample t-test
## 
## data:  df_long_acq_fix_count[df_long_acq_fix_count$iu_group == "1", "sticsa_total_centred"] and df_long_acq_fix_count[df_long_acq_fix_count$iu_group == "-1", "sticsa_total_centred"]
## t = 9.3255, df = 276, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   7.343247 11.273157
## sample estimates:
## mean of x mean of y 
##  4.754549 -4.553653
# p < .05 : sticsa is not independent of iu group 

# sticsa and stimulus
t_test_independence_sticsa_stimulus_acq_fix_count <- 
  t.test(
    df_long_acq_fix_count[df_long_acq_fix_count$stimulus == "1", "sticsa_total_centred"],
    df_long_acq_fix_count[df_long_acq_fix_count$stimulus == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_stimulus_acq_fix_count
## 
##  Two Sample t-test
## 
## data:  df_long_acq_fix_count[df_long_acq_fix_count$stimulus == "1", "sticsa_total_centred"] and df_long_acq_fix_count[df_long_acq_fix_count$stimulus == "-1", "sticsa_total_centred"]
## t = 0, df = 276, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.252832  2.252832
## sample estimates:
##                mean of x                mean of y 
## -0.000000000000002862893 -0.000000000000002862893
# p > .05 - sticsa is independent of stimulus

Extinction

# sticsa and iu group
t_test_independence_sticsa_iu_group_ext_fix_count <- 
  t.test(
    df_long_ext_fix_count[df_long_ext_fix_count$iu_group == "1", "sticsa_total_centred"],
    df_long_ext_fix_count[df_long_ext_fix_count$iu_group == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_iu_group_ext_fix_count
## 
##  Two Sample t-test
## 
## data:  df_long_ext_fix_count[df_long_ext_fix_count$iu_group == "1", "sticsa_total_centred"] and df_long_ext_fix_count[df_long_ext_fix_count$iu_group == "-1", "sticsa_total_centred"]
## t = 13.212, df = 554, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   7.924338 10.692067
## sample estimates:
## mean of x mean of y 
##  4.754549 -4.553653
# p < .05 : sticsa is not independent of iu group 

# sticsa and stimulus
t_test_independence_sticsa_stimulus_ext_fix_count <- 
  t.test(
    df_long_ext_fix_count[df_long_ext_fix_count$stimulus == "1", "sticsa_total_centred"],
    df_long_ext_fix_count[df_long_ext_fix_count$stimulus == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_stimulus_ext_fix_count
## 
##  Two Sample t-test
## 
## data:  df_long_ext_fix_count[df_long_ext_fix_count$stimulus == "1", "sticsa_total_centred"] and df_long_ext_fix_count[df_long_ext_fix_count$stimulus == "-1", "sticsa_total_centred"]
## t = 0, df = 554, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.586608  1.586608
## sample estimates:
##                mean of x                mean of y 
## -0.000000000000002862855 -0.000000000000002862855
# p > .05 - sticsa is independent of stimulus

# sticsa and time
t_test_independence_sticsa_time_ext_fix_count <- 
  t.test(
    df_long_ext_fix_count[df_long_ext_fix_count$time == "1", "sticsa_total_centred"],
    df_long_ext_fix_count[df_long_ext_fix_count$time == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_time_ext_fix_count
## 
##  Two Sample t-test
## 
## data:  df_long_ext_fix_count[df_long_ext_fix_count$time == "1", "sticsa_total_centred"] and df_long_ext_fix_count[df_long_ext_fix_count$time == "-1", "sticsa_total_centred"]
## t = 0, df = 554, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.586608  1.586608
## sample estimates:
##                mean of x                mean of y 
## -0.000000000000002862855 -0.000000000000002862855
# p > .05 - sticsa is independent of time

Fixation Duration (Log Transformed)

Acquisition

# sticsa and iu group
t_test_independence_sticsa_iu_group_acq_fix_duration_log <- 
  t.test(
    df_long_acq_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "1", "sticsa_total_centred"],
    df_long_acq_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_iu_group_acq_fix_duration_log
## 
##  Two Sample t-test
## 
## data:  df_long_acq_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "1", "sticsa_total_centred"] and df_long_acq_fix_duration_log[df_long_acq_fix_duration_log$iu_group == "-1", "sticsa_total_centred"]
## t = 9.3255, df = 276, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   7.343247 11.273157
## sample estimates:
## mean of x mean of y 
##  4.754549 -4.553653
# p < .05 : sticsa is not independent of iu group 

# sticsa and stimulus
t_test_independence_sticsa_stimulus_acq_fix_duration_log <- 
  t.test(
    df_long_acq_fix_duration_log[df_long_acq_fix_duration_log$stimulus == "1", "sticsa_total_centred"],
    df_long_acq_fix_duration_log[df_long_acq_fix_duration_log$stimulus == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_stimulus_acq_fix_duration_log
## 
##  Two Sample t-test
## 
## data:  df_long_acq_fix_duration_log[df_long_acq_fix_duration_log$stimulus == "1", "sticsa_total_centred"] and df_long_acq_fix_duration_log[df_long_acq_fix_duration_log$stimulus == "-1", "sticsa_total_centred"]
## t = 0, df = 276, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.252832  2.252832
## sample estimates:
##                mean of x                mean of y 
## -0.000000000000002862893 -0.000000000000002862893
# p > .05 - sticsa is independent of stimulus

Extinction

# sticsa and iu group
t_test_independence_sticsa_iu_group_ext_fix_duration_log <- 
  t.test(
    df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "1", "sticsa_total_centred"],
    df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_iu_group_ext_fix_duration_log
## 
##  Two Sample t-test
## 
## data:  df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "1", "sticsa_total_centred"] and df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$iu_group == "-1", "sticsa_total_centred"]
## t = 13.212, df = 554, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   7.924338 10.692067
## sample estimates:
## mean of x mean of y 
##  4.754549 -4.553653
# p < .05 : sticsa is not independent of iu group 

# sticsa and stimulus
t_test_independence_sticsa_stimulus_ext_fix_duration_log <- 
  t.test(
    df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$stimulus == "1", "sticsa_total_centred"],
    df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$stimulus == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_stimulus_ext_fix_duration_log
## 
##  Two Sample t-test
## 
## data:  df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$stimulus == "1", "sticsa_total_centred"] and df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$stimulus == "-1", "sticsa_total_centred"]
## t = 0, df = 554, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.586608  1.586608
## sample estimates:
##                mean of x                mean of y 
## -0.000000000000002862855 -0.000000000000002862855
# p > .05 - sticsa is independent of stimulus

# sticsa and time
t_test_independence_sticsa_time_ext_fix_duration <- 
  t.test(
    df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$time == "1", "sticsa_total_centred"],
    df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$time == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_time_ext_fix_duration
## 
##  Two Sample t-test
## 
## data:  df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$time == "1", "sticsa_total_centred"] and df_long_ext_fix_duration_log[df_long_ext_fix_duration_log$time == "-1", "sticsa_total_centred"]
## t = 0, df = 554, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.586608  1.586608
## sample estimates:
##                mean of x                mean of y 
## -0.000000000000002862855 -0.000000000000002862855
# p > .05 - sticsa is independent of time

Saccade Amplitude

Acquisition

# sticsa and iu group
t_test_independence_sticsa_iu_group_acq_sacc_amplitude <- 
  t.test(
    df_long_acq_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "1", "sticsa_total_centred"],
    df_long_acq_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_iu_group_acq_sacc_amplitude
## 
##  Two Sample t-test
## 
## data:  df_long_acq_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "1", "sticsa_total_centred"] and df_long_acq_sacc_amplitude[df_long_acq_sacc_amplitude$iu_group == "-1", "sticsa_total_centred"]
## t = 9.3255, df = 276, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   7.343247 11.273157
## sample estimates:
## mean of x mean of y 
##  4.754549 -4.553653
# p < .05 : sticsa is not independent of iu group 

# sticsa and stimulus
t_test_independence_sticsa_stimulus_acq_sacc_amplitude <- 
  t.test(
    df_long_acq_sacc_amplitude[df_long_acq_sacc_amplitude$stimulus == "1", "sticsa_total_centred"],
    df_long_acq_sacc_amplitude[df_long_acq_sacc_amplitude$stimulus == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_stimulus_acq_sacc_amplitude
## 
##  Two Sample t-test
## 
## data:  df_long_acq_sacc_amplitude[df_long_acq_sacc_amplitude$stimulus == "1", "sticsa_total_centred"] and df_long_acq_sacc_amplitude[df_long_acq_sacc_amplitude$stimulus == "-1", "sticsa_total_centred"]
## t = 0, df = 276, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.252832  2.252832
## sample estimates:
##                mean of x                mean of y 
## -0.000000000000002862893 -0.000000000000002862893
# p > .05 - sticsa is independent of stimulus

Extinction

# sticsa and iu group
t_test_independence_sticsa_iu_group_ext_sacc_amplitude <- 
  t.test(
    df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "1", "sticsa_total_centred"],
    df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_iu_group_ext_sacc_amplitude
## 
##  Two Sample t-test
## 
## data:  df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "1", "sticsa_total_centred"] and df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$iu_group == "-1", "sticsa_total_centred"]
## t = 13.212, df = 554, p-value < 0.00000000000000022
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   7.924338 10.692067
## sample estimates:
## mean of x mean of y 
##  4.754549 -4.553653
# p < .05 : sticsa is not independent of iu group 

# sticsa and stimulus
t_test_independence_sticsa_stimulus_ext_sacc_amplitude <- 
  t.test(
    df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$stimulus == "1", "sticsa_total_centred"],
    df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$stimulus == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_stimulus_ext_sacc_amplitude
## 
##  Two Sample t-test
## 
## data:  df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$stimulus == "1", "sticsa_total_centred"] and df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$stimulus == "-1", "sticsa_total_centred"]
## t = 0, df = 554, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.586608  1.586608
## sample estimates:
##                mean of x                mean of y 
## -0.000000000000002862855 -0.000000000000002862855
# p > .05 - sticsa is independent of stimulus

# sticsa and time
t_test_independence_sticsa_time_ext_sacc_amplitude <- 
  t.test(
    df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$time == "1", "sticsa_total_centred"],
    df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$time == "-1", "sticsa_total_centred"], 
    var.equal = TRUE
    )
t_test_independence_sticsa_time_ext_sacc_amplitude
## 
##  Two Sample t-test
## 
## data:  df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$time == "1", "sticsa_total_centred"] and df_long_ext_sacc_amplitude[df_long_ext_sacc_amplitude$time == "-1", "sticsa_total_centred"]
## t = 0, df = 554, p-value = 1
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.586608  1.586608
## sample estimates:
##                mean of x                mean of y 
## -0.000000000000002862855 -0.000000000000002862855
# p > .05 - sticsa is independent of time

Homogeneity of Regression Slopes

####### check homogeneity of regression slopes

  df_long_acq_fix_count %>%
    anova_test(fix_count ~ sticsa_total_centred + iu_group + stimulus + iu_group*stimulus + 
                 sticsa_total_centred*iu_group + sticsa_total_centred*stimulus +
                 sticsa_total_centred*iu_group*stimulus)
## Coefficient covariances computed by hccm()
## ANOVA Table (type II tests)
## 
##                                   Effect DFn DFd     F     p p<.05       ges
## 1                   sticsa_total_centred   1 270 0.114 0.736       0.0004210
## 2                               iu_group   1 270 6.146 0.014     * 0.0220000
## 3                               stimulus   1 270 0.957 0.329       0.0040000
## 4                      iu_group:stimulus   1 270 0.103 0.749       0.0003810
## 5          sticsa_total_centred:iu_group   1 270 3.336 0.069       0.0120000
## 6          sticsa_total_centred:stimulus   1 270 0.154 0.695       0.0005710
## 7 sticsa_total_centred:iu_group:stimulus   1 270 0.021 0.885       0.0000783
# as the interaction term sticsa_total_centred*iu_group was not sig (p = .067), there
  # is not an interaction between iu group and sticsa total centred
  
# as the interaction term sticsa_total_centred*stimulus was not sig (p = .787), there
  # is not an interaction between stimulus and sticsa total centred

  # There was homogeneity of regression slopes as the interaction terms between 
  # the covariate (STICSA total centred) and grouping variables (iu group
  # and stimulus), was not statistically significant, p > .05. 
  
# reporting example:
# There was homogeneity of regression slopes as the interaction term was not statistically significant, F(2, 39) = 0.13, p = 0.88.


  
######## check remaining assumptions
# to check remaining assumptions, first compute lm() version of the model
ancova_lm_model_acq_fix_count <- lm(fix_count ~ sticsa_total_centred + iu_group * stimulus, df_long_acq_fix_count)
model_metrics <- augment(ancova_lm_model_acq_fix_count) 
shapiro_test(model_metrics$.resid)
## # A tibble: 1 × 3
##   variable             statistic       p.value
##   <chr>                    <dbl>         <dbl>
## 1 model_metrics$.resid     0.940 0.00000000340
summary(ancova_lm_model_acq_fix_count)
## 
## Call:
## lm(formula = fix_count ~ sticsa_total_centred + iu_group * stimulus, 
##     data = df_long_acq_fix_count)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.4309 -2.6331 -0.4267  2.0982 15.4692 
## 
## Coefficients:
##                       Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)           6.706298   0.420831  15.936 <0.0000000000000002 ***
## sticsa_total_centred  0.008328   0.024721   0.337              0.7365    
## iu_group1             1.224693   0.623745   1.963              0.0506 .  
## stimulus1            -0.341613   0.573457  -0.596              0.5519    
## iu_group1:stimulus1  -0.120642   0.819886  -0.147              0.8831    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.417 on 273 degrees of freedom
## Multiple R-squared:  0.03625,    Adjusted R-squared:  0.02213 
## F-statistic: 2.567 on 4 and 273 DF,  p-value: 0.03851
# obtain the assumptions and save as variable 
q1_ancova_assumptions <- check_model(ancova_lm_model_acq_fix_count)

# check the assumptions
q1_ancova_assumptions
## Registered S3 methods overwritten by 'parameters':
##   method                           from      
##   as.double.parameters_kurtosis    datawizard
##   as.double.parameters_skewness    datawizard
##   as.double.parameters_smoothness  datawizard
##   as.numeric.parameters_kurtosis   datawizard
##   as.numeric.parameters_skewness   datawizard
##   as.numeric.parameters_smoothness datawizard
##   print.parameters_distribution    datawizard
##   print.parameters_kurtosis        datawizard
##   print.parameters_skewness        datawizard
##   summary.parameters_kurtosis      datawizard
##   summary.parameters_skewness      datawizard
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Loading required namespace: qqplotr

######## normality of residuals
# panel 2: most data points are roughly plotted along the line, but there do appear
# to be some possible deviations to assumption of normality 
# panel 3: residuals appear mostly normally distributed, 
# though there are some peaks on both tails
# assumption violated

###### homoscedasticity 
# panel 4: data-points roughly plotted along horizontal line
# assumption is met

####### outliers
# panel 6: no cases exceed Cook's value of 1, so no outliers by influence.
# assumption is met.

Outliers

## acquisition fix count
df_long_acq_fix_count %>%
  group_by(iu_group, stimulus) %>%
  identify_outliers(fix_count)
## # A tibble: 4 × 10
##   iu_group stimulus id    sticsa_total condition      fix_count sticsa_total_avg
##   <fct>    <fct>    <fct>        <dbl> <fct>              <dbl>            <dbl>
## 1 1        -1       086_1           68 acq_csm_fix_c…      18.3             40.5
## 2 1        -1       099_1           52 acq_csm_fix_c…      16               40.5
## 3 1        1        086_1           68 acq_csp_fix_c…      23.2             40.5
## 4 1        1        099_1           52 acq_csp_fix_c…      20.2             40.5
## # … with 3 more variables: sticsa_total_centred <dbl>, is.outlier <lgl>,
## #   is.extreme <lgl>
## extinction fix count
df_long_ext_fix_count %>%
  group_by(iu_group, stimulus, time) %>%
  identify_outliers(fix_count)
## # A tibble: 14 × 11
##    iu_group stimulus time  id    sticsa_total condition           fix_count
##    <fct>    <fct>    <fct> <fct>        <dbl> <fct>                   <dbl>
##  1 -1       -1       -1    122_1           37 l_ext_csm_fix_count      17.8
##  2 -1       -1       1     047_1           41 e_ext_csm_fix_count      17  
##  3 -1       -1       1     122_1           37 e_ext_csm_fix_count      21.5
##  4 -1       1        -1    122_1           37 l_ext_csp_fix_count      16  
##  5 -1       1        1     122_1           37 e_ext_csp_fix_count      20.5
##  6 -1       1        1     143_1           44 e_ext_csp_fix_count      19.5
##  7 1        -1       -1    033_1           54 l_ext_csm_fix_count      20  
##  8 1        -1       -1    065_1           33 l_ext_csm_fix_count      14.8
##  9 1        -1       -1    086_1           68 l_ext_csm_fix_count      19.2
## 10 1        -1       -1    099_1           52 l_ext_csm_fix_count      16  
## 11 1        -1       -1    113_1           31 l_ext_csm_fix_count      15  
## 12 1        1        -1    086_1           68 l_ext_csp_fix_count      22  
## 13 1        1        1     086_1           68 e_ext_csp_fix_count      19.2
## 14 1        1        1     113_1           31 e_ext_csp_fix_count      17.8
## # … with 4 more variables: sticsa_total_avg <dbl>, sticsa_total_centred <dbl>,
## #   is.outlier <lgl>, is.extreme <lgl>
# acquisition fix duration log
df_long_acq_fix_duration_log %>%
  group_by(iu_group, stimulus) %>%
  identify_outliers(fix_duration_log)
##  [1] iu_group             stimulus             id                  
##  [4] sticsa_total         condition            fix_duration_log    
##  [7] sticsa_total_avg     sticsa_total_centred is.outlier          
## [10] is.extreme          
## <0 rows> (or 0-length row.names)
## extinction fix count
df_long_ext_fix_duration_log %>%
  group_by(iu_group, stimulus, time) %>%
  identify_outliers(fix_duration_log)
## # A tibble: 6 × 11
##   iu_group stimulus time  id    sticsa_total condition          fix_duration_log
##   <fct>    <fct>    <fct> <fct>        <dbl> <fct>                         <dbl>
## 1 1        -1       -1    009_1           41 l_ext_csm_fix_dur…             8.00
## 2 1        -1       -1    010_1           43 l_ext_csm_fix_dur…             4.69
## 3 1        -1       -1    015_1           55 l_ext_csm_fix_dur…             8.37
## 4 1        -1       -1    044_1           36 l_ext_csm_fix_dur…             7.91
## 5 1        -1       1     044_1           36 e_ext_csm_fix_dur…             8.70
## 6 1        -1       1     113_1           31 e_ext_csm_fix_dur…             4.18
## # … with 4 more variables: sticsa_total_avg <dbl>, sticsa_total_centred <dbl>,
## #   is.outlier <lgl>, is.extreme <lgl>
## acquisition sacc amplitude
df_long_acq_sacc_amplitude %>%
  group_by(iu_group, stimulus) %>%
  identify_outliers(sacc_amplitude)
## # A tibble: 9 × 10
##   iu_group stimulus id    sticsa_total condition sacc_amplitude sticsa_total_avg
##   <fct>    <fct>    <fct>        <dbl> <fct>              <dbl>            <dbl>
## 1 -1       -1       016_1           26 acq_csm_…           7.37             40.5
## 2 -1       -1       026_1           55 acq_csm_…           6.87             40.5
## 3 -1       1        016_1           26 acq_csp_…           6.35             40.5
## 4 1        -1       017_1           33 acq_csm_…           7.81             40.5
## 5 1        -1       021_1           54 acq_csm_…           7.50             40.5
## 6 1        -1       022_1           50 acq_csm_…           8.57             40.5
## 7 1        1        009_1           41 acq_csp_…           7.47             40.5
## 8 1        1        043_1           39 acq_csp_…           6.88             40.5
## 9 1        1        044_1           36 acq_csp_…           8.15             40.5
## # … with 3 more variables: sticsa_total_centred <dbl>, is.outlier <lgl>,
## #   is.extreme <lgl>
## extinction sacc amplitude
df_long_ext_sacc_amplitude %>%
  group_by(iu_group, stimulus, time) %>%
  identify_outliers(sacc_amplitude)
## # A tibble: 17 × 11
##    iu_group stimulus time  id    sticsa_total condition           sacc_amplitude
##    <fct>    <fct>    <fct> <fct>        <dbl> <fct>                        <dbl>
##  1 -1       -1       -1    016_1           26 l_ext_csm_sacc_amp…          10.9 
##  2 -1       -1       -1    075_1           35 l_ext_csm_sacc_amp…           8.98
##  3 -1       -1       -1    078_1           42 l_ext_csm_sacc_amp…           8.03
##  4 -1       -1       -1    111_1           41 l_ext_csm_sacc_amp…           8.21
##  5 -1       -1       1     016_1           26 e_ext_csm_sacc_amp…           9.11
##  6 -1       1        -1    016_1           26 l_ext_csp_sacc_amp…          13.1 
##  7 -1       1        -1    075_1           35 l_ext_csp_sacc_amp…           7.84
##  8 -1       1        1     016_1           26 e_ext_csp_sacc_amp…           9.02
##  9 -1       1        1     051_1           28 e_ext_csp_sacc_amp…           7.40
## 10 -1       1        1     119_1           43 e_ext_csp_sacc_amp…           9.18
## 11 1        -1       -1    009_1           41 l_ext_csm_sacc_amp…           8.00
## 12 1        -1       -1    105_1           33 l_ext_csm_sacc_amp…           9.74
## 13 1        -1       1     105_1           33 e_ext_csm_sacc_amp…          11.4 
## 14 1        1        -1    009_1           41 l_ext_csp_sacc_amp…           9.62
## 15 1        1        -1    022_1           50 l_ext_csp_sacc_amp…           8.34
## 16 1        1        1     009_1           41 e_ext_csp_sacc_amp…           7.67
## 17 1        1        1     129_1           46 e_ext_csp_sacc_amp…           8.65
## # … with 4 more variables: sticsa_total_avg <dbl>, sticsa_total_centred <dbl>,
## #   is.outlier <lgl>, is.extreme <lgl>
### reporting:
# There were no extreme outliers.